diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright Tony Day (c) 2017 +Copyright Tony Day (c) 2016 All rights reserved. diff --git a/numhask-array.cabal b/numhask-array.cabal --- a/numhask-array.cabal +++ b/numhask-array.cabal @@ -1,5 +1,5 @@ name: numhask-array -version: 0.3.0.1 +version: 0.4.0.0 synopsis: n-dimensional arrays description: @@ -30,9 +30,7 @@ type: git location: - https://github.com/tonyday567/numhask - subdir: - numhask-array + https://github.com/tonyday567/numhask-array library hs-source-dirs: src @@ -51,17 +49,18 @@ base >=4.11 && <5 , adjunctions >=4.0 && <5 , deepseq >=1.4.2.0 && <2 - , dimensions >=1.0 && <1.1 , distributive >=0.4 && <0.7 + , numhask >= 0.3.1 && <0.4 , numhask-prelude >=0.3 && <0.4 , protolude >=0.2 && <0.3 - , singletons >=2.0 && <2.6 , vector >=0.10 && <0.13 + , hmatrix exposed-modules: NumHask.Array - NumHask.Array.Constraints - NumHask.Array.Example - NumHask.Shape + NumHask.Array.HMatrix + NumHask.Array.Fixed + NumHask.Array.Dynamic + NumHask.Array.Shape default-language: Haskell2010 test-suite test type: @@ -77,11 +76,10 @@ UnicodeSyntax build-depends: base >=4.11 && <5 + , adjunctions >=4.0 && <5 , doctest >=0.13 && <0.17 - , dimensions >=1.0 && <1.1 - , numhask-array >=0.3 && <0.4 + , numhask-array >=0.3 && <0.5 , numhask-prelude >=0.3 && <0.4 , numhask-hedgehog >=0.3 && <0.4 , hedgehog >=0.5 && <1.1 - , adjunctions >=4.0 && <5 default-language: Haskell2010 diff --git a/src/NumHask/Array.hs b/src/NumHask/Array.hs --- a/src/NumHask/Array.hs +++ b/src/NumHask/Array.hs @@ -1,707 +1,64 @@ -{-# LANGUAGE DataKinds #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE GADTs #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE NoImplicitPrelude #-} -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE PolyKinds #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE TypeApplications #-} -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE UndecidableInstances #-} -{-# OPTIONS_GHC -Wno-redundant-constraints #-} -{-# OPTIONS_GHC -fno-warn-missing-signatures #-} - -module NumHask.Array where - -import Data.Distributive (Distributive(..)) -import Data.Functor.Rep (Representable(..), liftR2, pureRep, fmapRep) -import Data.List ((!!)) -import GHC.Exts (IsList(..)) -import GHC.Show (Show(..)) -import NumHask.Error (impossible) -import NumHask.Array.Constraints - (Fold, HeadModule, TailModule, IsValidConcat, Concatenate, Transpose, Squeeze) -import NumHask.Prelude as P -import NumHask.Shape (HasShape(..)) -import Numeric.Dimensions as D -import qualified Data.Singletons.Prelude as S -import qualified Data.Vector as V - --- $setup --- >>> :set -XDataKinds --- >>> :set -XOverloadedLists --- >>> :set -XTypeFamilies --- >>> let a = [1..24] :: Array [] '[2,3,4] Int --- >>> let v = [1,2,3] :: Array [] '[3] Int +{-# OPTIONS_GHC -Wall #-} --- | an array polymorphic in container and shape +-- | Higher-kinded numbers that can be represented with an [Int] index into elements. -- --- >>> a --- [[[1, 2, 3, 4], --- [5, 6, 7, 8], --- [9, 10, 11, 12]], --- [[13, 14, 15, 16], --- [17, 18, 19, 20], --- [21, 22, 23, 24]]] -data family Array (c :: Type -> Type) (ds :: [k]) (a :: Type) - --- | instance where dimensions are known at compile time -newtype instance - Array c (ds :: [Nat]) t = - Array { _getContainer :: c t} - deriving (Functor, Foldable) - -instance NFData (Array c (ds :: [Nat]) t) where - rnf (Array c) = seq c () - -{- --- | instance of array where some of the dimensions are known at compile time --- it wraps an Array with some weird magic -data instance Array c (xds :: [XNat]) t = forall (ds :: [Nat]). - ( FixedDims xds ds - , Dimensions ds) => - SomeArray (Array c ds t) - --} - -instance (Dimensions r) => HasShape (Array c (r :: [Nat])) where - type Shape (Array c r) = [Int] - shape _ = fmap fromIntegral (listDims $ dims @Nat @r) - --- | an array with dimensions represented at the value level -newtype AnyArray c a = AnyArray ([Int], c a) - --- | convert an array with type-level shape to value-level shape -anyArray :: (Dimensions ds) => Array c (ds :: [Nat]) a -> AnyArray c a -anyArray arr@(Array c) = AnyArray (shape arr, c) - --- | a sweet class of container with attributes necessary to supply the set of operations here -class (Functor f) => Container f where - generate :: Int -> (Int -> a) -> f a - idx :: f a -> Int -> a - cslice :: Int -> Int -> f a -> f a - zipWith :: (a -> a -> a) -> f a -> f a -> f a - -- Chunks a container into a list of containers whose dimension are each i - chunkItUp :: [f a] -> Int -> f a -> [f a] - cfoldl' :: (b -> a -> b) -> b -> f a -> b - cfoldr :: (a -> b -> b) -> b -> f a -> b - cconcat :: [f a] -> f a - -instance Container V.Vector where - generate = V.generate - idx = V.unsafeIndex - cslice = V.unsafeSlice - zipWith = V.zipWith - chunkItUp acc i v = - if null v - then acc - else let (c, r) = V.splitAt i v - in chunkItUp (c : acc) i r - cfoldl' = V.foldl' - cfoldr = V.foldr - cconcat = V.concat - -instance Container [] where - generate n g = take n $ g <$> [0 ..] - idx = (!!) - cslice d t = take t . drop d - zipWith = P.zipWith - chunkItUp acc i v = - if null v - then acc - else let (c, r) = splitAt i v - in chunkItUp (c : acc) i r - cfoldl' = foldl' - cfoldr = foldr - cconcat = mconcat - -instance (Eq (c t), Dimensions ds) => Eq (Array c (ds :: [Nat]) t) where - (Array a) == (Array b) = a == b - -{- -dimList :: Dims ds -> [Int] -dimList U = [] -dimList (d :* ds) = dimList d ++ dimList ds -dimList (Dn _ :: Dim m) = [dimVal' @m] --- dimList (Dx (Dn _ :: Dim m)) = [dimVal' @m] +module NumHask.Array + ( -- * Shapes + -- + -- $shape + module NumHask.Array.Shape, --} + -- * Numbers with a fixed shape + -- + -- $fixed + module NumHask.Array.Fixed, + -- * Numbers with a dynamic shape + -- + -- $dynamic -{- -instance HasShape (Array c (xds :: [XNat])) where - type Shape (Array c xds) = [Int] - shape (SomeArray a) = shape a --} + -- * HMatrix API + -- + -- $hmatrix + ) where --- * shape helpers where dimensions ~ [Int] +import NumHask.Array.Shape +import NumHask.Array.Fixed --- | convert from n-dim shape index to a flat index +-- $shape -- --- >>> ind [2,3,4] [1,1,1] --- 17 -ind :: [Int] -> [Int] -> Int -ind ns xs = sum $ P.zipWith (*) xs (drop 1 $ scanr (*) 1 ns) +-- Using [`Int`] as the index for an array nicely represents the practical interests and constraints downstream of this high-level API: densely-packed numbers (reals or integrals), indexed and layered. --- | convert from a flat index to a shape index +-- $fixed -- --- >>> unind [2,3,4] 17 --- [1,1,1] -unind :: [Int] -> Int -> [Int] -unind ns x = - fst $ - foldr - (\a (acc, r) -> - let (d, m) = divMod r a - in (m : acc, d)) - ([], x) - ns - -instance forall r c. (Dimensions r, Container c) => - Data.Distributive.Distributive (Array c (r :: [Nat])) where - distribute f = Array $ generate (fromIntegral n) $ \i -> fmap (\(Array v) -> idx v i) f - where - n = totalDim $ dims @Nat @r - -instance forall r c. (Dimensions r, Container c) => - Representable (Array c (r :: [Nat])) where - type Rep (Array c r) = [Int] - tabulate f = Array $ generate (fromIntegral $ product ns) (f . unind (fmap fromIntegral ns)) - where - ns = listDims $ dims @Nat @r - index (Array xs) rs = xs `idx` ind (fmap fromIntegral ns) rs - where - ns = listDims $ dims @Nat @r - --- | from flat list -instance - ( Item (Array c r a) ~ Item (c a) - , Dimensions r - , Additive a - , IsList (c a) - ) => - IsList (Array c (r :: [Nat]) a) where - type Item (Array c r a) = a - fromList l = Array $ fromList $ take n $ l ++ repeat zero - where - n = fromIntegral $ totalDim (dims @_ @r) - toList (Array v) = GHC.Exts.toList v - -instance (Show a, Show (Item (c a)), Container c, IsList (c a)) => Show (AnyArray c a) where - show aa@(AnyArray (l,_)) = go (length l) aa - where - go n aa'@(AnyArray (l', c')) = - case length l' of - 0 -> "[]" - 1 -> "[" ++ intercalate ", " (GHC.Show.show <$> GHC.Exts.toList c') ++ "]" - x -> - "[" ++ - intercalate - (",\n" ++ replicate (n - x + 1) ' ') - (go n <$> flatten1 aa') ++ - "]" - --- | convert the top layer of a SomeArray to a [SomeArray] -flatten1 :: (Container c) => AnyArray c a -> [AnyArray c a] -flatten1 (AnyArray (rep, v)) = - (\s -> AnyArray (drop 1 rep, cslice (s * l) l v)) <$> ss - where - (n, l) = - case rep of - [] -> (0, 1) - x:r -> (x, product r) - ss = take n [0 ..] - -instance (Show a, Show (Item (c a)), IsList (c a), Container c, Dimensions ds) - => Show (Array c (ds :: [Nat]) a) where - show = GHC.Show.show . anyArray - -type Vector c n = Array c '[ n] - -type Matrix c m n = Array c '[ m, n] - -{- -instance - ( IsList (c a) - , Item (c a) ~ a - , Container c - , KnownNat n - , Unital (Sum (Vector c n a)) - , QC.Arbitrary a - , Additive a) => - QC.Arbitrary (Vector c n a) where - arbitrary = QC.frequency [(1, pure zero), (9, fromList <$> QC.vector n)] - where - n = fromInteger $ P.natVal (Proxy :: Proxy n) - -instance - ( IsList (c a) - , Item (c a) ~ a - , Additive (Matrix c m n a) - , Container c - , KnownNat m - , KnownNat n - , QC.Arbitrary a - , Additive a) => - QC.Arbitrary (Matrix c m n a) where - arbitrary = QC.frequency [(1, pure zero), (9, fromList <$> QC.vector (m * n))] - where - n = fromInteger $ P.natVal (Proxy :: Proxy n) - m = fromInteger $ P.natVal (Proxy :: Proxy m) --} +-- A design principle of modern haskell (synonymous with ghc usage) is to push computation into compilation. `Representable` is an iconic example. `Shape` must ultimately be known at compile-time, but if it is, we can dispense with a large classes of runtime check which slow pipelines. This tends to create computation where shape is abstracted and we fold algorithms in and out of structures. --- ** Operations --- | outer product +-- $dynamic -- --- todo: reconcile with numhask version +-- | In many situations, where shape is being tracked or otherwise only known at runtime, a clear module arrangement is: -- --- >>> v NumHask.Array.>< v --- [[1, 2, 3], --- [2, 4, 6], --- [3, 6, 9]] -(><) :: forall c (r :: [Nat]) (s :: [Nat]) a. - ( Container c - , CommutativeRing a - , Dimensions r - , Dimensions s - , Dimensions ((D.++) r s)) - => Array c r a - -> Array c s a - -> Array c ((D.++) r s) a -(><) m n = tabulate (\i -> index m (take dimm i) * index n (drop dimm i)) - where - dimm = length (shape m) +-- >>> import NumHask.Array.Shape +-- >>> import qualified NumHask.Array.Fixed as F +-- >>> import qualified NumHask.Array.Dynamic as D --- | matrix multiplication +-- $hmatrix -- --- >>> let a = [1, 2, 3, 4] :: Array [] '[2, 2] Int --- >>> let b = [5, 6, 7, 8] :: Array [] '[2, 2] Int --- >>> a --- [[1, 2], --- [3, 4]] +-- | hmatrix remains the speed king within haskell, and numhask-array is no exception. If speed matters then: -- --- >>> b --- [[5, 6], --- [7, 8]] +-- >>> import NumHask.Array.Shape +-- >>> import qualified NumHask.Array.HMatrix as H -- --- >>> mmult a b --- [[19, 22], --- [43, 50]] +-- will deliver the NumHask API over a HMatrix representation, and then the current fastest way to multiply matices in boiler-plate haskell. -- -mmult :: forall c m n k a. - ( Hilbert (Vector c k a) - , Dimensions '[ m, k] - , Dimensions '[ k, n] - , Dimensions '[ m, n] - , Container c - ) - => Matrix c (m :: Nat) (k :: Nat) a - -> Matrix c k n a - -> Matrix c m n a -mmult x y = tabulate go - where - go [i, j] = unsafeRow i x <.> unsafeCol j y - go _ = impossible "mmult only typechecks on arrays" +-- FIXME: fill out Fixed API --- | extract the row of a matrix -row :: forall c i a m n. - ( Dimensions '[ m, n] - , Container c - , KnownNat i - , ((S.<) i m) ~ 'True - ) - => Proxy i - -> Matrix c m n a - -> Vector c n a -row i_ = unsafeRow i - where - i = (fromIntegral . S.fromSing . S.singByProxy) i_ -rank2Shape - :: Dimensions '[ m, n] - => Matrix c (m :: Nat) (n :: Nat) a - -> (Int, Int) -rank2Shape t = - case shape t of - [m, n] -> (m, n) - _ -> impossible "only typechecks for matricies" -unsafeRow :: forall c a m n. - ( Container c - , Dimensions '[ m, n]) - => Int - -> Matrix c (m :: Nat) (n :: Nat) a - -> Vector c n a -unsafeRow i t@(Array a) = Array $ cslice (i * n) n a - where - (_, n) = rank2Shape t -unsafeCol :: - forall c a m n. (Container c, Dimensions '[ m, n]) - => Int - -> Matrix c (m :: Nat) (n :: Nat) a - -> Vector c m a -unsafeCol j t@(Array a) = Array $ generate m (\x -> a `idx` (j + x * n)) - where - (m, n) = rank2Shape t --- | extract the column of a matrix -col :: forall c j a m n. - ( Dimensions '[ m, n] - , Container c - , KnownNat j - , ((S.<) j n) ~ 'True - ) - => Proxy j - -> Matrix c m n a - -> Vector c m a -col j_ = unsafeCol j - where - j = (fromIntegral . S.fromSing . S.singByProxy) j_ --- | --- --- >>> unsafeIndex a [0,2,1] --- 10 -unsafeIndex :: (Container c, Dimensions r) => Array c (r :: [Nat]) a -> [Int] -> a -unsafeIndex t@(Array a) i = a `idx` ind (shape t) i --- | --- --- >>> unsafeSlice [[0,1],[2],[1,2]] a :: Array [] '[2,1,2] Int --- [[[10, 11]], --- [[22, 23]]] -unsafeSlice :: - (Container c, IsList (c a), Item (c a) ~ a, Dimensions r, Dimensions r0) - => [[Int]] - -> Array c (r :: [Nat]) a - -> Array c (r0 :: [Nat]) a -unsafeSlice s t = Array (fromList [unsafeIndex t i | i <- sequence s]) --- | --- --- todo: an ambiguous type variable has snuck in here somewhere --- --- > slice (Proxy :: Proxy '[ '[0,1],'[2],'[1,2]]) a --- [[[10, 11]], --- [[22, 23]]] -{- -todo: - • Expected kind ‘[[Nat]]’, but ‘s’ has kind ‘[Nat]’ - • In the first argument of ‘Slice’, namely ‘s’ - In the first argument of ‘Array’, namely ‘(Slice s)’ - In the type signature: - slice :: forall c s r a. - (Container c, - Dimensions s, - Dimensions r, - And (ZipWith AllLTSym0 s r) ~ 'True) => - Proxy s -> Array c r a -> Array (Slice s) c a --} -{- -slice :: - forall c s r a. (Container c, Dimensions s, Dimensions r, S.And (S.ZipWith AllLTSym0 s r) ~ 'True) - => Proxy s - -> Array c r a - -> Array (Slice s) c a --} -slice s_ = unsafeSlice s - where - s = ((fmap . fmap) fromInteger . S.fromSing . S.singByProxy) s_ --- | --- --- >>> foldAlong (Proxy :: Proxy 1) (\_ -> ([0..3] :: Array [] '[4] Int)) a --- [[0, 1, 2, 3], --- [0, 1, 2, 3]] --- --- todo: resolution of a primitive and a scalar eg --- Expected type: Array '[10] Int -> Array '[] Int --- Actual type: Array '[10] (Array '[] Int) -> Array '[] Int --- -foldAlong :: - forall c s vw uvw uw w a. - ( Container c - , KnownNat s - , Dimensions uvw - , uw ~ (Fold s uvw) - , w ~ (S.Drop 1 vw) - , vw ~ (TailModule s uvw) - ) - => Proxy s - -> (Array c vw a -> Array c w a) - -> Array c uvw a - -> Array c uw a -foldAlong s_ f a@(Array v) = - Array $ - cconcat - (cfoldl' - (\xs x -> - let (Array vx) = f (Array x) - in vx : xs) - [] - md) - where - s = (fromIntegral . S.fromSing . S.singByProxy) s_ - md = chunkItUp [] (product $ drop s $ shape a) v - --- | --- --- todo: No instance for (Container (Array [] '[]) error --- --- > mapAlong (Proxy :: Proxy 0) (\x -> NumHask.Array.zipWith (*) x x) a --- [[[1, 4, 9, 16], --- [25, 36, 49, 64], --- [81, 100, 121, 144]], --- [[169, 196, 225, 256], --- [289, 324, 361, 400], --- [441, 484, 529, 576]]] --- -mapAlong :: - forall c s uvw vw a. - (Container c, KnownNat s, Dimensions uvw, vw ~ (HeadModule s uvw)) - => Proxy s - -> (Array c vw a -> Array c vw a) - -> Array c uvw a - -> Array c uvw a -mapAlong s_ f a@(Array v) = - Array $ - cconcat - (cfoldl' - (\xs x -> - let (Array vx) = f (Array x) - in vx : xs) - [] - md) - where - s = (fromIntegral . S.fromSing . S.singByProxy) s_ - md = chunkItUp [] (product $ drop s $ shape a) v - --- | --- --- >>> concatenate (Proxy :: Proxy 2) a a --- [[[1, 2, 3, 4, 1, 2, 3, 4], --- [5, 6, 7, 8, 5, 6, 7, 8], --- [9, 10, 11, 12, 9, 10, 11, 12]], --- [[13, 14, 15, 16, 13, 14, 15, 16], --- [17, 18, 19, 20, 17, 18, 19, 20], --- [21, 22, 23, 24, 21, 22, 23, 24]]] --- -concatenate :: - forall c s r t a. - ( Container c - , S.SingI s - , Dimensions r - , Dimensions t - , (IsValidConcat s t r) ~ 'True - ) - => Proxy s - -> Array c r a - -> Array c t a - -> Array c (Concatenate s t r) a -concatenate s_ r@(Array vr) t@(Array vt) = - Array . cconcat $ (concat . reverse . P.transpose) [rm, tm] - where - s = (fromIntegral . S.fromSing . S.singByProxy) s_ - rm = chunkItUp [] (product $ drop s $ shape t) vt - tm = chunkItUp [] (product $ drop s $ shape r) vr - --- | --- --- >>> NumHask.Array.transpose a --- [[[1, 2], --- [3, 4], --- [5, 6]], --- [[7, 8], --- [9, 10], --- [11, 12]], --- [[13, 14], --- [15, 16], --- [17, 18]], --- [[19, 20], --- [21, 22], --- [23, 24]]] --- -transpose :: - forall c s t a. (t ~ Transpose s, Container c, Dimensions s, Dimensions t) - => Array c (s :: [Nat]) a - -> Array c (t :: [Nat]) a -transpose (Array x) = Array x - --- | --- --- >>> let a = [1..24] :: Array [] '[2,1,3,4,1] Int --- >>> a --- [[[[[1], --- [2], --- [3], --- [4]], --- [[5], --- [6], --- [7], --- [8]], --- [[9], --- [10], --- [11], --- [12]]]], --- [[[[13], --- [14], --- [15], --- [16]], --- [[17], --- [18], --- [19], --- [20]], --- [[21], --- [22], --- [23], --- [24]]]]] --- >>> squeeze a --- [[[1, 2, 3, 4], --- [5, 6, 7, 8], --- [9, 10, 11, 12]], --- [[13, 14, 15, 16], --- [17, 18, 19, 20], --- [21, 22, 23, 24]]] --- -squeeze :: - forall c s t a. (t ~ Squeeze s) - => Array c s a - -> Array c t a -squeeze (Array x) = Array x - -instance (Dimensions r, Container c, Additive a) => - Additive (Array c (r :: [Nat]) a) where - a + b = liftR2 (+) a b - zero = pureRep zero - -instance (Dimensions r, Container c, Subtractive a) => - Subtractive (Array c (r :: [Nat]) a) where - negate = fmapRep negate - -instance (Dimensions r, Container c, Multiplicative a) => - Multiplicative (Array c (r :: [Nat]) a) where - a * b = liftR2 (*) a b - - one = pureRep one - -instance (Dimensions r, Container c, Divisive a) => - Divisive (Array c (r :: [Nat]) a) where - recip = fmapRep recip - -instance (Dimensions r, Container c, Multiplicative a, Additive a) => - P.Distributive (Array c (r :: [Nat]) a) - -instance (Dimensions r, Container c, IntegralDomain a) => IntegralDomain (Array c (r :: [Nat]) a) - -instance (Dimensions r, Container c, Field a) => Field (Array c (r :: [Nat]) a) - -instance (Dimensions r, Container c, ExpField a) => ExpField (Array c (r :: [Nat]) a) where - exp = fmapRep exp - log = fmapRep log - -instance (Foldable (Array c r), Dimensions r, Container c, UpperBoundedField a) => - UpperBoundedField (Array c (r :: [Nat]) a) where - isNaN = foldl' (||) False . fmapRep isNaN - -instance (Foldable (Array c r), Dimensions r, Container c, LowerBoundedField a) => - LowerBoundedField (Array c (r :: [Nat]) a) - -instance (Dimensions r, Container c, Multiplicative a, Signed a) - => Signed (Array c (r :: [Nat]) a) where - sign = fmapRep sign - abs = fmapRep abs - -instance (Functor (Array c r), Foldable (Array c r), Additive (Array c r a), Normed a a, ExpField a) => - Normed (Array c (r :: [Nat]) a) a where - normL1 r = foldr (+) zero $ normL1 <$> r - normL2 r = sqrt $ foldr (+) zero $ (** (one + one)) <$> r - -instance (Eq (c a), Foldable (Array c r), Dimensions r, Container c, Epsilon a) => - Epsilon (Array c (r :: [Nat]) a) where - epsilon = tabulate (const epsilon) - nearZero f = and (fmapRep nearZero f) - aboutEqual a b = and (liftR2 aboutEqual a b) - -instance (Foldable (Array c r), Dimensions r, Container c, ExpField a, Subtractive a, Normed a a) => - Metric (Array c (r :: [Nat]) a) a where - distanceL1 a b = normL1 (a - b) - distanceL2 a b = normL2 (a - b) - -instance (Dimensions r, Container c, Integral a) => Integral (Array c (r :: [Nat]) a) where - divMod a b = (d, m) - where - x = liftR2 divMod a b - d = fmap fst x - m = fmap snd x - quotRem a b = (q, r) - where - x = liftR2 quotRem a b - q = fmap fst x - r = fmap snd x - -type instance Actor (Array c r a) = a - -instance (Dimensions r, Container c, Multiplicative a) => - HadamardMultiplication (Array c (r :: [Nat])) a where - (.*.) = liftR2 (*) - -instance (Dimensions r, Container c, Divisive a) => - HadamardDivision (Array c (r :: [Nat])) a where - (./.) = liftR2 (/) - -instance (Dimensions r, Container c, Additive a) => - AdditiveAction (Array c (r::[Nat]) a) where - (.+) r s = fmap (s +) r - (+.) s = fmap (s +) - -instance (Dimensions r, Container c, Subtractive a) => - SubtractiveAction (Array c (r::[Nat]) a) where - (.-) r s = fmap (\x -> x - s) r - (-.) s = fmap (\x -> x - s) - -instance (Dimensions r, Container c, Multiplicative a) => - MultiplicativeAction (Array c (r :: [Nat]) a) where - (.*) r s = fmap (* s) r - (*.) s = fmap (s *) - -instance (Dimensions r, Container c, Divisive a) => - DivisiveAction (Array c (r::[Nat]) a) where - (./) r s = fmap (/ s) r - (/.) s = fmap (/ s) - -instance forall a c r. (Actor (Array c r a) ~ a, Foldable (Array c r), P.Distributive a, CommutativeRing a, Semiring a, Dimensions r, Container c) => - Hilbert (Array c (r :: [Nat]) a) where - a <.> b = sum $ liftR2 (*) a b - -instance - ( Foldable (Array c r) - , Dimensions r - , Container c - , CommutativeRing a - , Multiplicative a - ) => - TensorProduct (Array c (r :: [Nat]) a) where - (><) m n = tabulate (\i -> index m i *. n) - timesleft v m = tabulate (\i -> v <.> index m i) - timesright m v = tabulate (\i -> v <.> index m i) - -instance (Eq (c a), Container c, Dimensions r, JoinSemiLattice a) => JoinSemiLattice (Array c (r :: [Nat]) a) where - (\/) = liftR2 (\/) - -instance (Eq (c a), Container c, Dimensions r, MeetSemiLattice a) => MeetSemiLattice (Array c (r :: [Nat]) a) where - (/\) = liftR2 (/\) - -instance (Eq (c a), Container c, Dimensions r, BoundedJoinSemiLattice a) => BoundedJoinSemiLattice (Array c (r :: [Nat]) a) where - bottom = pureRep bottom - -instance (Eq (c a), Container c, Dimensions r, BoundedMeetSemiLattice a) => BoundedMeetSemiLattice (Array c (r :: [Nat]) a) where - top = pureRep top - -singleton :: (Dimensions r, Container c) => a -> Array c (r :: [Nat]) a -singleton a = tabulate (const a) diff --git a/src/NumHask/Array/Constraints.hs b/src/NumHask/Array/Constraints.hs deleted file mode 100644 --- a/src/NumHask/Array/Constraints.hs +++ /dev/null @@ -1,83 +0,0 @@ -{-# LANGUAGE CPP #-} -{-# LANGUAGE ConstraintKinds #-} -{-# LANGUAGE DataKinds #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE GADTs #-} -{-# LANGUAGE NoImplicitPrelude #-} -{-# LANGUAGE PolyKinds #-} -{-# LANGUAGE RankNTypes #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE TypeInType #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE UndecidableInstances #-} -{-# OPTIONS_GHC -Wall #-} -{-# OPTIONS_GHC -fno-warn-deprecations #-} -{-# OPTIONS_GHC -fno-warn-orphans #-} - -module NumHask.Array.Constraints - ( IsValidConcat - , Squeeze - , Concatenate - , IsValidTranspose - , Fold - , FoldAlong - , TailModule - , HeadModule - , Transpose - ) where - -import Data.Singletons.Prelude hiding (Max) -import Data.Singletons.Prelude.List hiding (Transpose) -import Data.Singletons.Prelude.Tuple (Fst, Snd) -import Data.Singletons.TypeLits (Nat) -import qualified Protolude as P - -type family DropDim d a :: [b] where - DropDim 0 xs = Drop 1 xs - DropDim d xs = Take (d - 1) (Fst (SplitAt d xs)) ++ Snd (SplitAt d xs) - -type family IsValidConcat i (a :: [Nat]) (b :: [Nat]) :: P.Bool where - IsValidConcat _ '[] _ = 'P.False - IsValidConcat _ _ '[] = 'P.False - IsValidConcat i a b = - And (ZipWith (==@#@$) (DropDim i a) (DropDim i b)) - -type family Squeeze (a :: [Nat]) where - Squeeze '[] = '[] - Squeeze a = Filter ((/=@#@$$) 1) a - -type family IsValidTranspose (p :: [Nat]) (a :: [Nat]) :: P.Bool where - IsValidTranspose p a = - (Minimum p >= 0) && (Minimum a >= 0) && (Sum a == Sum p) && Length p == Length a - -type family Transpose a where - Transpose a = Reverse a - -type family AddDimension (d :: Nat) t :: [Nat] where - AddDimension d t = Insert d t - -type family Concatenate i (a :: [Nat]) (b :: [Nat]) :: [Nat] where - Concatenate i a b = - Take i (Fst (SplitAt (i + 1) a)) ++ - ('[ Head (Drop i a) + Head (Drop i b)]) ++ - Snd (SplitAt (i + 1) b) - --- | Reduces axis i in shape s. Maintains singlton dimension -type family FoldAlong i (s :: [Nat]) where - FoldAlong _ '[] = '[] - FoldAlong d xs = Take d (Fst (SplitAt (d + 1) xs)) ++ '[ 1] ++ Snd (SplitAt (d + 1) xs) - --- | Reduces axis i in shape s. Does not maintain singlton dimension. -type family Fold i (s :: [Nat]) where - Fold _ '[] = '[] - Fold d xs = Take d (Fst (SplitAt (d + 1) xs)) ++ Snd (SplitAt (d + 1) xs) - -type family TailModule i (s :: [Nat]) where - TailModule _ '[] = '[] - TailModule d xs = (Snd (SplitAt d xs)) - -type family HeadModule i (s :: [Nat]) where - HeadModule _ '[] = '[] - HeadModule d xs = (Fst (SplitAt d xs)) diff --git a/src/NumHask/Array/Dynamic.hs b/src/NumHask/Array/Dynamic.hs new file mode 100644 --- /dev/null +++ b/src/NumHask/Array/Dynamic.hs @@ -0,0 +1,548 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE DeriveAnyClass #-} +{-# LANGUAGE DeriveGeneric #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE NoStarIsType #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# OPTIONS_GHC -Wno-redundant-constraints #-} +{-# OPTIONS_GHC -Wno-incomplete-uni-patterns #-} + +-- | Arrays with a dynamic shape +module NumHask.Array.Dynamic + ( -- * Dynamic Arrays + -- + -- $array + Array(..), + fromFlatList, + + -- * Operators + -- + -- $operators + reshape, + transpose, + diag, + selects, + selectsExcept, + folds, + extracts, + joins, + maps, + concatenate, + insert, + append, + reorder, + expand, + contract, + dot, + slice, + squeeze, + singleton, + ident, + -- * Scalar + -- + -- $scalar + fromScalar, + toScalar, + -- * Matrix + -- + -- $matrix + col, + row, + mmult, + )where + +import qualified Data.Vector as V +import GHC.Show (Show (..)) +import NumHask.Array.Shape +import NumHask.Prelude as P hiding (product, transpose) +import Data.List ((!!)) + +-- $setup +-- >>> :set -XDataKinds +-- >>> :set -XOverloadedLists +-- >>> :set -XTypeFamilies +-- >>> :set -XFlexibleContexts +-- >>> let s = fromFlatList [] [1] :: Array Int +-- >>> let a = fromFlatList [2,3,4] [1..24] :: Array Int +-- >>> let v = fromFlatList [3] [1,2,3] :: Array Int +-- >>> let m = fromFlatList [3,4] [0..11] :: Array Int + +-- | a multidimensional array with a value-level shape +-- +-- >>> a +-- [[[1, 2, 3, 4], +-- [5, 6, 7, 8], +-- [9, 10, 11, 12]], +-- [[13, 14, 15, 16], +-- [17, 18, 19, 20], +-- [21, 22, 23, 24]]] +data Array a + = Array {shape :: [Int], unArray :: V.Vector a} + deriving (Eq, Ord, NFData, Generic) + +instance Functor Array where + fmap f (Array s a) = Array s (V.map f a) + +instance Foldable Array where + foldr x a (Array _ v) = V.foldr x a v + +instance Traversable Array where + traverse f (Array s v) = + fromFlatList s <$> traverse f (toList v) + +instance (Show a) => Show (Array a) where + show a@(Array l _) = go (length l) a + where + go n a'@(Array l' m) = + case length l' of + 0 -> maybe (throw (NumHaskException "empty scalar")) GHC.Show.show (head m) + 1 -> "[" ++ intercalate ", " (GHC.Show.show <$> V.toList m) ++ "]" + x -> + "[" + ++ intercalate + (",\n" ++ replicate (n - x + 1) ' ') + (go n <$> toFlatList (extracts [0] a')) + ++ "]" + +-- * conversions + +-- | convert from a list +-- +-- >>> fromFlatList [2,3,4] [1..24] == a +-- True +fromFlatList :: [Int] -> [a] -> Array a +fromFlatList ds l = Array ds $ V.fromList $ take (size ds) l + +-- | convert to a flat list. +-- +-- >>> toFlatList a == [1..24] +-- True +toFlatList :: Array a -> [a] +toFlatList (Array _ v) = V.toList v + +-- | extract an element at index /i/ +-- +-- >>> index a [1,2,3] +-- 24 +index :: () => Array a -> [Int] -> a +index (Array s v) i = V.unsafeIndex v (flatten s i) + +-- | tabulate an array with a generating function +-- +-- >>> tabulate [2,3,4] ((1+) . flatten [2,3,4]) == a +-- True +tabulate :: () => [Int] -> ([Int] -> a) -> Array a +tabulate ds f = Array ds . V.generate (size ds) $ (f . shapen ds) + +-- | reshape an array (with the same number of elements) +-- +-- >>> reshape [4,3,2] a +-- [[[1, 2], +-- [3, 4], +-- [5, 6]], +-- [[7, 8], +-- [9, 10], +-- [11, 12]], +-- [[13, 14], +-- [15, 16], +-- [17, 18]], +-- [[19, 20], +-- [21, 22], +-- [23, 24]]] +reshape :: + [Int] -> + Array a -> + Array a +reshape s a = tabulate s (index a . shapen (shape a) . flatten s) + +-- | reverse indices eg transposes the element /Aijk/ to /Akji/ +-- +-- >>> index (transpose a) [1,0,0] == index a [0,0,1] +-- True +transpose :: Array a -> Array a +transpose a = tabulate (reverse $ shape a) (index a . reverse) + +-- | +-- +-- >>> ident [3,2] +-- [[1, 0], +-- [0, 1], +-- [0, 0]] +ident :: (Num a) => [Int] -> Array a +ident ds = tabulate ds (bool 0 1 . isDiag) + where + isDiag [] = True + isDiag [_] = True + isDiag [x, y] = x == y + isDiag (x : y : xs) = x == y && isDiag (y : xs) + +-- | extract the diagonal +-- +-- >>> diag (ident [3,2]) +-- [1, 1] +diag :: + Array a -> + Array a +diag a = tabulate [NumHask.Array.Shape.minimum (shape a)] go + where + go [] = throw (NumHaskException "Rank Underflow") + go (s' : _) = index a (replicate (rank (shape a)) s') + +-- | +-- >>> singleton [3,2] one +-- [[1, 1], +-- [1, 1], +-- [1, 1]] +singleton :: [Int] -> a -> Array a +singleton ds a = tabulate ds (const a) + +-- | /selects ds ps a/ select from /a/, elements along /ds/ dimensions at positions /ps/ +-- +-- >>> let s = selects [0,1] [1,1] a +-- >>> s +-- [17, 18, 19, 20] +selects :: + [Int] -> + [Int] -> + Array a -> + Array a +selects ds i a = tabulate (dropIndexes (shape a) ds) go + where + go s = index a (addIndexes s ds i) + +-- | select an index /except/ along dimensions +-- +-- >>> let s = selectsExcept [2] [1,1] a +-- >>> s +-- [17, 18, 19, 20] +selectsExcept :: + [Int] -> + [Int] -> + Array a -> + Array a +selectsExcept ds i a = selects (exclude (rank (shape a)) ds) i a + +-- | fold along specified dimensions +-- +-- >>> folds sum [1] a +-- [68, 100, 132] +folds :: + (Array a -> b) -> + [Int] -> + Array a -> + Array b +folds f ds a = tabulate (takeIndexes (shape a) ds) go + where + go s = f (selects ds s a) + +-- | extracts dimensions to an outer layer +-- +-- >>> let e = extracts [1,2] a +-- >>> shape <$> extracts [0] a +-- [[3,4], [3,4]] +extracts :: + [Int] -> + Array a -> + Array (Array a) +extracts ds a = tabulate (takeIndexes (shape a) ds) go + where + go s = selects ds s a + +-- | extracts /except/ dimensions to an outer layer +-- +-- >>> let e = extractsExcept [1,2] a +-- >>> shape <$> extracts [0] a +-- [[3,4], [3,4]] +extractsExcept :: + [Int] -> + Array a -> + Array (Array a) +extractsExcept ds a = extracts (exclude (rank (shape a)) ds) a + +-- | join inner and outer dimension layers +-- +-- >>> let e = extracts [1,0] a +-- >>> let j = joins [1,0] e +-- >>> a == j +-- True +joins :: + [Int] -> + Array (Array a) -> + Array a +joins ds a = tabulate (addIndexes si ds so) go + where + go s = index (index a (takeIndexes s ds)) (dropIndexes s ds) + so = shape a + si = shape (index a (replicate (rank so) 0)) + +-- | maps along specified dimensions +-- +-- >>> shape $ maps (transpose) [1] a +-- [4,3,2] +maps :: + (Array a -> Array b) -> + [Int] -> + Array a -> + Array b +maps f ds a = joins ds (fmap f (extracts ds a)) + +-- | concatenate along a dimension +-- +-- >>> shape $ concatenate 1 a a +-- [2,6,4] +concatenate :: + Int -> + Array a -> + Array a -> + Array a +concatenate d a0 a1 = tabulate (concatenate' d (shape a0) (shape a1)) go + where + go s = + bool + (index a0 s) + ( index + a1 + ( addIndex + (dropIndex s d) + d + ((s !! d) - (ds0 !! d)) + ) + ) + ((s !! d) >= (ds0 !! d)) + ds0 = shape a0 + +-- | /insert d i/ insert along the dimension /d/ at position /i/ +-- +-- >>> insert 2 0 a (fromFlatList [2,3] [100..105]) +-- [[[100, 1, 2, 3, 4], +-- [101, 5, 6, 7, 8], +-- [102, 9, 10, 11, 12]], +-- [[103, 13, 14, 15, 16], +-- [104, 17, 18, 19, 20], +-- [105, 21, 22, 23, 24]]] +insert :: + Int -> + Int -> + Array a -> + Array a -> + Array a +insert d i a b = tabulate (incAt d (shape a)) go + where + go s + | s !! d == i = index b (dropIndex s d) + | s !! d < i = index a s + | otherwise = index a (decAt d s) + +-- | insert along a dimension at the end +-- +-- >>> append 2 a (fromFlatList [2,3] [100..105]) +-- [[[1, 2, 3, 4, 100], +-- [5, 6, 7, 8, 101], +-- [9, 10, 11, 12, 102]], +-- [[13, 14, 15, 16, 103], +-- [17, 18, 19, 20, 104], +-- [21, 22, 23, 24, 105]]] +append :: + Int -> + Array a -> + Array a -> + Array a +append d a b = insert d (dimension (shape a) d) a b + +-- | change the order of dimensions +-- +-- >>> let r = reorder [2,0,1] a +-- >>> r +-- [[[1, 5, 9], +-- [13, 17, 21]], +-- [[2, 6, 10], +-- [14, 18, 22]], +-- [[3, 7, 11], +-- [15, 19, 23]], +-- [[4, 8, 12], +-- [16, 20, 24]]] +reorder :: + [Int] -> + Array a -> + Array a +reorder ds a = tabulate (reorder' (shape a) ds) go + where + go s = index a (addIndexes [] ds s) + +-- | product two arrays using the supplied binary function +-- If the function is multiply, and the arrays are tensors, +-- then this can be interpreted as a tensor product. +-- +-- https://en.wikipedia.org/wiki/Tensor_product +-- +-- The concept of a tensor product is a dense crossroad, and a complete treatment is elsewhere. To quote: +-- ... the tensor product can be extended to other categories of mathematical objects in addition to vector spaces, such as to matrices, tensors, algebras, topological vector spaces, and modules. In each such case the tensor product is characterized by a similar universal property: it is the freest bilinear operation. The general concept of a "tensor product" is captured by monoidal categories; that is, the class of all things that have a tensor product is a monoidal category. +-- +-- >>> expand (*) v v +-- [[1, 2, 3], +-- [2, 4, 6], +-- [3, 6, 9]] +expand :: + (a -> b -> c) -> + Array a -> + Array b -> + Array c +expand f a b = tabulate ((++) (shape a) (shape b)) (\i -> f (index a (take r i)) (index b (drop r i))) + where + r = rank (shape a) + +-- | contract an array by applying the supplied (folding) function on diagonal elements of the dimensions. +-- +-- This generalises a tensor contraction by allowing the number of contracting diagonals to be other than 2, and allowing another binary other than addition +-- +-- >>> let b = fromFlatList [2,3] [1..6] :: Array Int +-- >>> contract sum [1,2] (expand (*) b (transpose b)) +-- [[14, 32], +-- [32, 77]] +contract :: + (Array a -> b) -> + [Int] -> + Array a -> + Array b +contract f xs a = f . diag <$> extractsExcept xs a + +-- | a generalisation of a dot operation, which is a multiplicative expansion of two arrays and sum contraction along the middle two dimensions. +-- +-- dot sum (*) on two matrices is known as matrix multiplication +-- +-- >>> let b = fromFlatList [2,3] [1..6] :: Array Int +-- >>> dot sum (*) b (transpose b) +-- [[14, 32], +-- [32, 77]] +-- +-- dot sum (*) on two vectors is known as the inner product +-- +-- >>> let v = fromFlatList [3] [1..3] :: Array Int +-- >>> dot sum (*) v v +-- 14 +-- +-- dot sum (*) m v on a matrix and a vector is matrix-vector multiplication +-- Note that an `Array Int` with shape [3] is neither a row vector nor column vector. `dot` is not turning the vector into a matrix and then using matrix multiplication. +-- +-- >>> dot sum (*) v b +-- [9, 12, 15] +-- +-- >>> dot sum (*) b v +-- [14, 32] +dot :: + (Array c -> d) -> + (a -> b -> c) -> + Array a -> + Array b -> + Array d +dot f g a b = contract f [rank sa - 1, rank sa] (expand g a b) + where + sa = shape a + +-- | select elements along every dimension +-- +-- >>> let s = slice [[0,1],[0,2],[1,2]] a +-- >>> s +-- [[[2, 3], +-- [10, 11]], +-- [[14, 15], +-- [22, 23]]] +-- +slice :: + [[Int]] -> + Array a -> + Array a +slice pss a = tabulate (ranks pss) go + where + go s = index a (zipWith (!!) pss s) + +-- | remove singleton dimensions +-- +-- >>> let a' = fromFlatList [2,1,3,4,1] [1..24] :: Array Int +-- >>> shape $ squeeze a' +-- [2,3,4] +squeeze :: + Array a -> + Array a +squeeze (Array s x) = Array (squeeze' s) x + +-- * scalar specializations + +-- | +-- +-- An /Array/ with shape [] despite being a Scalar is nevertheless a one-element vector under the hood. + +-- | unwrapping scalars is probably a performance bottleneck +-- +-- >>> let s = fromFlatList [] [3] :: Array Int +-- >>> fromScalar s +-- 3 +fromScalar :: Array a -> a +fromScalar a = index a ([] :: [Int]) + +-- | convert a number to a scalar +-- +-- >>> :t toScalar 2 +-- toScalar 2 :: Num a => Array a +toScalar :: a -> Array a +toScalar a = fromFlatList [] [a] + +-- * matrix specializations +-- | extract specialised to a matrix +-- +-- >>> row 1 m +-- [4, 5, 6, 7] +row :: Int -> Array a -> Array a +row i (Array s a) = Array [n] $ V.slice (i * n) n a + where + [_,n] = s + +-- | extract specialised to a matrix +-- +-- >>> col 1 m +-- [1, 5, 9] +col :: Int -> Array a -> Array a +col i (Array s a) = Array [m] $ V.generate m (\x -> V.unsafeIndex a (i + x * n)) + where + [m,n] = s + +-- | matrix multiplication +-- +-- This is dot sum (*) specialised to matrices +-- +-- >>> let a = fromFlatList [2,2] [1, 2, 3, 4] :: Array Int +-- >>> let b = fromFlatList [2,2] [5, 6, 7, 8] :: Array Int +-- >>> a +-- [[1, 2], +-- [3, 4]] +-- +-- >>> b +-- [[5, 6], +-- [7, 8]] +-- +-- >>> mmult a b +-- [[19, 22], +-- [43, 50]] +mmult :: + (Ring a) => + Array a -> + Array a -> + Array a +mmult (Array sx x) (Array sy y) = tabulate [m,n] go + where + go [] = throw (NumHaskException "Needs two dimensions") + go [_] = throw (NumHaskException "Needs two dimensions") + go (i : j : _) = sum $ V.zipWith (*) (V.slice (fromIntegral i * k) k x) (V.generate k (\x' -> y V.! (fromIntegral j + x' * n))) + [m,k] = sx + [_,n] = sy +{-# INLINE mmult #-} + + diff --git a/src/NumHask/Array/Example.hs b/src/NumHask/Array/Example.hs deleted file mode 100644 --- a/src/NumHask/Array/Example.hs +++ /dev/null @@ -1,250 +0,0 @@ -{-# LANGUAGE OverloadedLists #-} -{-# LANGUAGE DataKinds #-} -{-# LANGUAGE NoImplicitPrelude #-} -{-# OPTIONS_GHC -Wall #-} -{-# OPTIONS_GHC -fno-warn-unused-imports #-} - --- | Experimental api following https://pechersky.github.io/haskell-numpy-docs/quickstart.basics.html -module NumHask.Array.Example - ( - -- * The Basics - -- $setup - - -- ** An Example - -- $anExample - - -- ** Array Creation - -- $arrayCreation - - -- ** Printing Arrays - -- $printingArrays - - -- ** Universal Functions - -- $universalFunctions - - -- ** Indexing, Slicing and Iterating - -- $indexingSlicingIterating - - -- * Shape Manipulation - -- $shapeManipulation - - -- * Fancy indexing and index tricks - -- $fancyIndexing - - -- * Linear Algebra - -- $linearAlgebra - - -- * Tricks and Tips - -- $tricksTips - ) -where - -import NumHask.Shape -import NumHask.Prelude as P - --- $setup --- >>> :set -XDataKinds --- >>> :set -XOverloadedLists --- >>> :set -XNoImplicitPrelude --- >>> :set -XFlexibleContexts --- >>> import NumHask.Array as A --- >>> import GHC.Exts (fromList) - --- $anExample --- construction can be lazy; and zero pads --- --- >>> let z = [] :: Array [] '[2] Int --- >>> z --- [0, 0] --- >>> let a = [0..] :: Array [] '[3,5] Int --- >>> a --- [[0, 1, 2, 3, 4], --- [5, 6, 7, 8, 9], --- [10, 11, 12, 13, 14]] --- >>> shape a --- [3,5] --- >>> length (shape a) -- dimension --- 2 --- >>> :t a --- a :: Array [] '[3, 5] Int --- >>> import qualified Data.Vector as V --- >>> let v = V.fromList [6,7,8] --- >>> :t v --- v :: Num a => V.Vector a --- >>> let b = Array v --- >>> :t b --- b :: Num t => Array V.Vector ds t --- >>> b :: Array V.Vector '[3] Int --- [6, 7, 8] - --- $arrayCreation --- --- >>> -- fixed size arrays are fully shape specified at the type level --- >>> let a = [2, 3, 4] :: Array [] '[3] Int --- >>> a --- [2, 3, 4] --- >>> [1.2, 3.5, 5.1] :: Array [] '[3] Double --- [1.2, 3.5, 5.1] --- --- >>> -- lists of lists is not a thing, and need to be flattened --- >>> let ls = [[1.0,2.0,3.0],[4.0,5.0,6.0]] :: [[Double]] --- >>> fromList (concat ls) :: Array [] '[2,3] Double --- [[1.0, 2.0, 3.0], --- [4.0, 5.0, 6.0]] --- --- >>> fromList ((\x -> (fromIntegral x) :+ zero) <$> [1,2,3,4]) :: Array [] '[2,2] (Complex Double) --- [[1.0 :+ 0.0, 2.0 :+ 0.0], --- [3.0 :+ 0.0, 4.0 :+ 0.0]] --- >>> let z = [] :: Array [] '[3,4] Int --- >>> z --- [[0, 0, 0, 0], --- [0, 0, 0, 0], --- [0, 0, 0, 0]] --- >>> let o = A.singleton one :: Array [] '[2,3,4] Int --- >>> o --- [[[1, 1, 1, 1], --- [1, 1, 1, 1], --- [1, 1, 1, 1]], --- [[1, 1, 1, 1], --- [1, 1, 1, 1], --- [1, 1, 1, 1]]] --- >>> let empt = A.singleton nan :: Array [] '[2,3] Double --- >>> empt --- [[NaN, NaN, NaN], --- [NaN, NaN, NaN]] --- --- >>> [10,15 .. 30] :: Array [] '[4] Int --- [10, 15, 20, 25] --- >>> [0, 0.3.. 2] :: Array [] '[7] Double --- [0.0, 0.3, 0.6, 0.8999999999999999, 1.2, 1.5, 1.7999999999999998] --- --- > todo: fix NumHask.Range grid --- > fromList (grid OuterPos (Range 0 2) 8) :: Array [] '[9] Double --- > [0.0, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0] --- > let x = fromList (grid OuterPos (Range 0 (2*pi)) 100) :: Array [] '[101] Double --- > let f = fmap sin x --- - --- $printingArrays --- >>> show ([0..] :: Array [] '[6] Int) :: Text --- "[0, 1, 2, 3, 4, 5]" --- >>> [0..] :: Array [] '[2,3] Int --- [[0, 1, 2], --- [3, 4, 5]] --- >>> [0..] :: Array [] '[1,2,3,1] Int --- [[[[0], --- [1], --- [2]], --- [[3], --- [4], --- [5]]]] --- --- > todo: implement display --- > import Formatting --- > display (left 7 . fixed 2) ", " (fromList $ fromIntegral <$> [0..] :: Array [] '[100,100] Double) --- > [[ 0.00, 1.00, 2.00 .. 98.00 99.00], --- > [ 100.00, 101.00, 102.00 .. 198.00 199.00], --- > [ 200.00, 201.00, 202.00 .. 298.00 299.00], --- > .. --- > [9800.00, 9801.00, 9802.00 .. 9898.00 9899.00], --- > [9900.00, 9901.00, 9902.00 .. 9998.00 9999.00]] --- - --- $basicOperation --- --- >>> let a = [20,30,40,50] :: Array [] '[4] Double --- >>> let b = [0..] :: Array [] '[4] Double --- >>> let c = a - b --- >>> c --- [20.0, 29.0, 38.0, 47.0] --- >>> -- todo: resolve potential to polymorph number literals eg b**2 --- >>> b ** (one+one) --- [0.0, 1.0, 4.0, 9.000000000000002] --- >>> 10 *. (sin <$> a) --- [9.129452507276277, -9.880316240928618, 7.451131604793488, -2.6237485370392877] --- >>> (<35) <$> a --- [True, True, False, False] --- --- >>> let a = [0..] :: Array [] '[2,3] Int --- >>> let b = [0..] :: Array [] '[3,2] Int --- >>> a .*. a --- [[0, 1, 4], --- [9, 16, 25]] --- >>> -- todo: resolve NumHask.Array and NumHask.Matrix operation names --- >>> a `NumHask.Array.mmult` b --- [[10, 13], --- [28, 40]] --- >>> a .* 2 --- [[0, 2, 4], --- [6, 8, 10]] --- --- >>> -- random example skipped --- --- > -- todo: awaiting grid fix --- > let a = singleton one :: Array [] '[3] Double --- > let b = fromList (grid OuterPos (Range 0 pi) 2) :: Array [] '[3] Double --- > let c = a + b --- > let d = exp . ((zero:+one) *.) $ (:+zero) <$> c --- > d --- [0.5403023058681398 :+ 0.8414709848078965, (-0.8414709848078965) :+ 0.5403023058681398, (-0.5403023058681399) :+ (-0.8414709848078964)] --- > :t d --- d :: Array [] '[3] (Complex Double) --- --- >>> -- folding --- >>> let a = [0..] :: Array [] '[2,5] Int --- >>> sum a --- 45 --- >>> minimum a --- 0 --- >>> maximum a --- 9 --- >>> -- todo: scanAlong --- - --- $universalFunctions --- --- >>> let a = [0..] :: Array [] '[3] Double --- >>> exp <$> a --- [1.0, 2.718281828459045, 7.38905609893065] --- >>> sqrt <$> a --- [0.0, 1.0, 1.4142135623730951] --- - --- $indexingSlicingIterating --- --- >>> let a = (\x -> x*x*x) <$> [0..] :: Array [] '[10] Int --- >>> index a [2] --- 8 --- >>> let s = (\i -> index a [i]) <$> [2..5] :: Array [] '[4] Int --- >>> s --- [8, 27, 64, 125] --- >>> :t s --- s :: Array [] '[4] Int --- >>> -- replace every second number with -1000 --- >>> let a' = (tabulate (\[i] -> if i `mod` 2 == 0 then -1000 else (index a [i]))) :: Array [] '[10] Int --- >>> a' --- [-1000, 1, -1000, 27, -1000, 125, -1000, 343, -1000, 729] --- --- > -- todo: reverse fix --- > let a'' = (let (n:_) = shape a in tabulate (\[i] -> index a [n-i])) :: Array [] '[4] Int --- > a'' --- > [729, -1000, 343, -1000, 125, -1000, 27, -1000, 1, -1000] --- --- > -- todo: slicing api - --- $shapeManipulation --- --- > -- todo: - --- $fancyIndexing --- --- > -- todo: - --- $linearAlgebra --- --- > -- todo: - --- $tricksTips --- --- > -- todo: - diff --git a/src/NumHask/Array/Fixed.hs b/src/NumHask/Array/Fixed.hs new file mode 100644 --- /dev/null +++ b/src/NumHask/Array/Fixed.hs @@ -0,0 +1,911 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE DeriveFoldable #-} +{-# LANGUAGE DeriveFunctor #-} +{-# LANGUAGE DeriveGeneric #-} +{-# LANGUAGE DeriveTraversable #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE NoStarIsType #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeApplications #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# OPTIONS_GHC -Wno-redundant-constraints #-} + +-- | Arrays with a fixed shape. +module NumHask.Array.Fixed + ( -- * Fixed-sixed arrays + -- + -- $array + Array (..), + with, + shape, + toDynamic, + + -- * Operators + -- + -- $operators + reshape, + transpose, + diag, + selects, + selectsExcept, + folds, + extracts, + joins, + maps, + concatenate, + insert, + append, + reorder, + expand, + contract, + dot, + slice, + squeeze, + ident, + singleton, + + -- * Scalar + -- + -- $scalar + Scalar, + fromScalar, + toScalar, + + -- * Vector + -- + -- $scalar + Vector, + + -- * Matrix + -- + -- $matrix + Matrix, + col, + row, + safeCol, + safeRow, + mmult, + ) +where + +import Data.Distributive (Distributive (..)) +import Data.Functor.Rep +import Data.List ((!!)) +import qualified Data.Vector as V +import GHC.Exts (IsList (..)) +import GHC.Show (Show (..)) +import GHC.TypeLits +import qualified NumHask.Array.Dynamic as D +import NumHask.Array.Shape +import NumHask.Prelude as P hiding (identity, outer, transpose) + +-- $setup +-- >>> :set -XDataKinds +-- >>> :set -XOverloadedLists +-- >>> :set -XTypeFamilies +-- >>> :set -XFlexibleContexts +-- >>> let s = [1] :: Array ('[] :: [Nat]) Int +-- >>> let a = [1..24] :: Array '[2,3,4] Int +-- >>> let v = [1,2,3] :: Array '[3] Int +-- >>> let m = [0..11] :: Array '[3,4] Int + +-- | a multidimensional array with a type-level shape +-- +-- >>> a +-- [[[1, 2, 3, 4], +-- [5, 6, 7, 8], +-- [9, 10, 11, 12]], +-- [[13, 14, 15, 16], +-- [17, 18, 19, 20], +-- [21, 22, 23, 24]]] +-- >>> [1,2,3] :: Array '[2,2] Int +-- [[*** Exception: NumHaskException {errorMessage = "shape mismatch"} +newtype Array s a = Array {unArray :: V.Vector a} deriving (Eq, Ord, NFData, Functor, Foldable, Generic, Traversable) + +instance (HasShape s, Show a) => Show (Array s a) where + show a = GHC.Show.show (toDynamic a) + +instance + ( HasShape s + ) => + Data.Distributive.Distributive (Array s) + where + distribute = distributeRep + {-# INLINE distribute #-} + +instance + forall s. + ( HasShape s + ) => + Representable (Array s) + where + + type Rep (Array s) = [Int] + + tabulate f = + Array . V.generate (size s) $ (f . shapen s) + where + s = shapeVal $ toShape @s + {-# INLINE tabulate #-} + + index (Array v) i = V.unsafeIndex v (flatten s i) + where + s = shapeVal (toShape @s) + {-# INLINE index #-} + +-- * NumHask heirarchy + +instance + ( Additive a, + HasShape s + ) => + Additive (Array s a) + where + + (+) = liftR2 (+) + + zero = pureRep zero + +instance + ( Subtractive a, + HasShape s + ) => + Subtractive (Array s a) + where + negate = fmapRep negate + +type instance Actor (Array s a) = a + +instance + ( Multiplicative a, + HasShape s + ) => + HadamardMultiplication (Array s) a + where + (.*.) = liftR2 (*) + +instance + ( Divisive a, + HasShape s + ) => + HadamardDivision (Array s) a + where + (./.) = liftR2 (/) + +instance + (HasShape s, Multiplicative a) => + MultiplicativeAction (Array s a) + where + + (.*) r s = fmap (* s) r + {-# INLINE (.*) #-} + + (*.) s = fmap (s *) + {-# INLINE (*.) #-} + +instance (HasShape s, JoinSemiLattice a) => JoinSemiLattice (Array s a) where + (\/) = liftR2 (\/) + +instance (HasShape s, MeetSemiLattice a) => MeetSemiLattice (Array s a) where + (/\) = liftR2 (/\) + +instance (HasShape s, Subtractive a, Epsilon a) => Epsilon (Array s a) where + + epsilon = singleton epsilon + + nearZero (Array a) = all nearZero a + +-- | from flat list +instance + ( HasShape s + ) => + IsList (Array s a) + where + + type Item (Array s a) = a + + fromList l = + bool + (throw (NumHaskException "shape mismatch")) + (Array $ V.fromList l) + ((length l == 1 && null ds) || (length l == size ds)) + where + ds = shapeVal (toShape @s) + + toList (Array v) = V.toList v + +-- * shape + +-- | get shape of an Array as a value +-- +-- >>> shape a +-- [2,3,4] +shape :: forall a s. (HasShape s) => Array s a -> [Int] +shape _ = shapeVal $ toShape @s +{-# INLINE shape #-} + +-- | convert to a dynamic array with shape at the value level. +toDynamic :: (HasShape s) => Array s a -> D.Array a +toDynamic a = D.fromFlatList (shape a) (P.toList a) + +-- | use a dynamic array in a fixed context +-- +-- >>> with (D.fromFlatList [2,3,4] [1..24]) (selects (Proxy :: Proxy '[0,1]) [1,1] :: Array '[2,3,4] Int -> Array '[4] Int) +-- [17, 18, 19, 20] +with :: + forall a r s. + (HasShape s) => + D.Array a -> + (Array s a -> r) -> + r +with (D.Array _ v) f = f (Array v) + +-- * operations + +-- | reshape an array (with the same number of elements) +-- +-- >>> reshape a :: Array '[4,3,2] Int +-- [[[1, 2], +-- [3, 4], +-- [5, 6]], +-- [[7, 8], +-- [9, 10], +-- [11, 12]], +-- [[13, 14], +-- [15, 16], +-- [17, 18]], +-- [[19, 20], +-- [21, 22], +-- [23, 24]]] +reshape :: + forall a s s'. + ( Size s ~ Size s', + HasShape s, + HasShape s' + ) => + Array s a -> + Array s' a +reshape a = tabulate (index a . shapen s . flatten s') + where + s = shapeVal (toShape @s) + s' = shapeVal (toShape @s') + +-- | reverse indices eg transposes the element /Aijk/ to /Akji/ +-- +-- >>> index (transpose a) [1,0,0] == index a [0,0,1] +-- True +transpose :: forall a s. (HasShape s, HasShape (Reverse s)) => Array s a -> Array (Reverse s) a +transpose a = tabulate (index a . reverse) + +-- | +-- +-- >>> ident :: Array '[3,2] Int +-- [[1, 0], +-- [0, 1], +-- [0, 0]] +ident :: forall a s. (HasShape s, Additive a, Multiplicative a) => Array s a +ident = tabulate (bool zero one . isDiag) + where + isDiag [] = True + isDiag [_] = True + isDiag [x, y] = x == y + isDiag (x : y : xs) = x == y && isDiag (y : xs) + +-- | +-- >>> diag (ident :: Array '[3,2] Int) +-- [1, 1] +diag :: + forall a s. + ( HasShape s, + HasShape '[Minimum s] + ) => + Array s a -> + Array '[Minimum s] a +diag a = tabulate go + where + go [] = throw (NumHaskException "Rank Underflow") + go (s' : _) = index a (replicate (length ds) s') + ds = shapeVal (toShape @s) + +-- | +-- >>> singleton one :: Array '[3,2] Int +-- [[1, 1], +-- [1, 1], +-- [1, 1]] +singleton :: (HasShape s) => a -> Array s a +singleton a = tabulate (const a) + +-- | /selects ds ps a/ select from /a/ elements /ds/ dimensions at positions /ps/ +-- +-- >>> let s = selects (Proxy :: Proxy '[0,1]) [1,1] a +-- >>> :t s +-- s :: Array '[4] Int +-- +-- >>> s +-- [17, 18, 19, 20] +selects :: + forall ds s s' a. + ( HasShape s, + HasShape ds, + HasShape s', + s' ~ DropIndexes s ds + ) => + Proxy ds -> + [Int] -> + Array s a -> + Array s' a +selects _ i a = tabulate go + where + go s = index a (addIndexes s ds i) + ds = shapeVal (toShape @ds) + +-- | select an index /except/ along dimensions +-- +-- >>> let s = selectsExcept (Proxy :: Proxy '[2]) [1,1] a +-- >>> :t s +-- s :: Array '[4] Int +-- +-- >>> s +-- [17, 18, 19, 20] +selectsExcept :: + forall ds s s' a. + ( HasShape s, + HasShape ds, + HasShape s', + s' ~ TakeIndexes s ds + ) => + Proxy ds -> + [Int] -> + Array s a -> + Array s' a +selectsExcept _ i a = tabulate go + where + go s = index a (addIndexes i ds s) + ds = shapeVal (toShape @ds) + +-- | fold along specified dimensions +-- +-- >>> folds sum (Proxy :: Proxy '[1]) a +-- [68, 100, 132] +folds :: + forall ds st si so a b. + ( HasShape st, + HasShape ds, + HasShape si, + HasShape so, + si ~ DropIndexes st ds, + so ~ TakeIndexes st ds + ) => + (Array si a -> b) -> + Proxy ds -> + Array st a -> + Array so b +folds f d a = tabulate go + where + go s = f (selects d s a) + +-- | extracts dimensions to an outer layer +-- +-- >>> let e = extracts (Proxy :: Proxy '[1,2]) a +-- >>> :t e +-- e :: Array '[3, 4] (Array '[2] Int) +extracts :: + forall ds st si so a. + ( HasShape st, + HasShape ds, + HasShape si, + HasShape so, + si ~ DropIndexes st ds, + so ~ TakeIndexes st ds + ) => + Proxy ds -> + Array st a -> + Array so (Array si a) +extracts d a = tabulate go + where + go s = selects d s a + +-- | extracts /except/ dimensions to an outer layer +-- +-- >>> let e = extractsExcept (Proxy :: Proxy '[1,2]) a +-- >>> :t e +-- e :: Array '[2] (Array '[3, 4] Int) +extractsExcept :: + forall ds st si so a. + ( HasShape st, + HasShape ds, + HasShape si, + HasShape so, + so ~ DropIndexes st ds, + si ~ TakeIndexes st ds + ) => + Proxy ds -> + Array st a -> + Array so (Array si a) +extractsExcept d a = tabulate go + where + go s = selectsExcept d s a + +-- | join inner and outer dimension layers +-- +-- >>> let e = extracts (Proxy :: Proxy '[1,0]) a +-- +-- >>> :t e +-- e :: Array '[3, 2] (Array '[4] Int) +-- +-- >>> let j = joins (Proxy :: Proxy '[1,0]) e +-- +-- >>> :t j +-- j :: Array '[2, 3, 4] Int +-- +-- >>> a == j +-- True +joins :: + forall ds si st so a. + ( HasShape st, + HasShape ds, + st ~ AddIndexes si ds so, + HasShape si, + HasShape so + ) => + Proxy ds -> + Array so (Array si a) -> + Array st a +joins _ a = tabulate go + where + go s = index (index a (takeIndexes s ds)) (dropIndexes s ds) + ds = shapeVal (toShape @ds) + +-- | maps along specified dimensions +-- +-- >>> :t maps (transpose) (Proxy :: Proxy '[1]) a +-- maps (transpose) (Proxy :: Proxy '[1]) a :: Array '[4, 3, 2] Int +maps :: + forall ds st st' si si' so a b. + ( HasShape st, + HasShape st', + HasShape ds, + HasShape si, + HasShape si', + HasShape so, + si ~ DropIndexes st ds, + so ~ TakeIndexes st ds, + st' ~ AddIndexes si' ds so, + st ~ AddIndexes si ds so + ) => + (Array si a -> Array si' b) -> + Proxy ds -> + Array st a -> + Array st' b +maps f d a = joins d (fmapRep f (extracts d a)) + +-- | concatenate along a dimension +-- +-- >>> :t concatenate (Proxy :: Proxy 1) a a +-- concatenate (Proxy :: Proxy 1) a a :: Array '[2, 6, 4] Int +concatenate :: + forall a s0 s1 d s. + ( CheckConcatenate d s0 s1 s, + Concatenate d s0 s1 ~ s, + HasShape s0, + HasShape s1, + HasShape s, + KnownNat d + ) => + Proxy d -> + Array s0 a -> + Array s1 a -> + Array s a +concatenate _ s0 s1 = tabulate go + where + go s = + bool + (index s0 s) + ( index + s1 + ( addIndex + (dropIndex s d) + d + ((s !! d) - (ds0 !! d)) + ) + ) + ((s !! d) >= (ds0 !! d)) + ds0 = shapeVal (toShape @s0) + d = fromIntegral $ natVal @d Proxy + +-- | /insert (Proxy :: Proxy d) (Proxy :: Proxy i)/ insert along the dimension /d/ at position /i/ +-- +-- >>> insert (Proxy :: Proxy 2) (Proxy :: Proxy 0) a ([100..105]) +-- [[[100, 1, 2, 3, 4], +-- [101, 5, 6, 7, 8], +-- [102, 9, 10, 11, 12]], +-- [[103, 13, 14, 15, 16], +-- [104, 17, 18, 19, 20], +-- [105, 21, 22, 23, 24]]] +insert :: + forall a s s' d i. + ( DropIndex s d ~ s', + CheckInsert d i s, + KnownNat i, + KnownNat d, + HasShape s, + HasShape s', + HasShape (Insert d s) + ) => + Proxy d -> + Proxy i -> + Array s a -> + Array s' a -> + Array (Insert d s) a +insert _ _ a b = tabulate go + where + go s + | s !! d == i = index b (dropIndex s d) + | s !! d < i = index a s + | otherwise = index a (decAt d s) + d = fromIntegral $ natVal @d Proxy + i = fromIntegral $ natVal @i Proxy + +-- | insert along a dimension at the end +-- +-- >>> :t append (Proxy :: Proxy 0) a +-- append (Proxy :: Proxy 0) a +-- :: Array '[3, 4] Int -> Array '[3, 3, 4] Int +append :: + forall a d s s'. + ( DropIndex s d ~ s', + CheckInsert d (Dimension s d - 1) s, + KnownNat (Dimension s d - 1), + KnownNat d, + HasShape s, + HasShape s', + HasShape (Insert d s) + ) => + Proxy d -> + Array s a -> + Array s' a -> + Array (Insert d s) a +append d = insert d (Proxy :: Proxy (Dimension s d - 1)) + +-- | change the order of dimensions +-- +-- >>> let r = reorder (Proxy :: Proxy '[2,0,1]) a +-- >>> :t r +-- r :: Array '[4, 2, 3] Int +reorder :: + forall a ds s. + ( HasShape ds, + HasShape s, + HasShape (Reorder s ds), + CheckReorder ds s + ) => + Proxy ds -> + Array s a -> + Array (Reorder s ds) a +reorder _ a = tabulate go + where + go s = index a (addIndexes [] ds s) + ds = shapeVal (toShape @ds) + +-- | product two arrays using the supplied binary function +-- If the function is multiply, and the arrays are tensors, +-- then this can be interpreted as a tensor product. +-- +-- https://en.wikipedia.org/wiki/Tensor_product +-- +-- The concept of a tensor product is a dense crossroad, and a complete treatment is elsewhere. To quote: +-- ... the tensor product can be extended to other categories of mathematical objects in addition to vector spaces, such as to matrices, tensors, algebras, topological vector spaces, and modules. In each such case the tensor product is characterized by a similar universal property: it is the freest bilinear operation. The general concept of a "tensor product" is captured by monoidal categories; that is, the class of all things that have a tensor product is a monoidal category. +-- +-- >>> expand (*) v v +-- [[1, 2, 3], +-- [2, 4, 6], +-- [3, 6, 9]] +expand :: + forall s s' a b c. + ( HasShape s, + HasShape s', + HasShape ((++) s s') + ) => + (a -> b -> c) -> + Array s a -> + Array s' b -> + Array ((++) s s') c +expand f a b = tabulate (\i -> f (index a (take r i)) (index b (drop r i))) + where + r = rank (shape a) + +-- | contract an array by applying the supplied (folding) function on diagonal elements of the dimensions. +-- +-- This generalises a tensor contraction by allowing the number of contracting diagonals to be other than 2, and allowing another binary other than addition +-- +-- >>> let b = [1..6] :: Array '[2,3] Int +-- >>> contract sum (Proxy :: Proxy '[1,2]) (expand (*) b (transpose b)) +-- [[14, 32], +-- [32, 77]] +contract :: + forall a b s ss s' ds. + ( KnownNat (Minimum (TakeIndexes s ds)), + HasShape (TakeIndexes s ds), + HasShape s, + HasShape ds, + HasShape ss, + HasShape s', + s' ~ DropIndexes s ds, + ss ~ '[Minimum (TakeIndexes s ds)] + ) => + (Array ss a -> b) -> + Proxy ds -> + Array s a -> + Array s' b +contract f xs a = f . diag <$> extractsExcept xs a + +-- | a generalisation of a dot operation, which is a multiplicative expansion of two arrays and sum contraction along the middle two dimensions. +-- +-- dot sum (*) on two matrices is known as matrix multiplication +-- +-- >>> let b = [1..6] :: Array '[2,3] Int +-- >>> dot sum (*) b (transpose b) +-- [[14, 32], +-- [32, 77]] +-- +-- dot sum (*) on two vectors is known as the inner product +-- +-- >>> let v = [1..3] :: Array '[3] Int +-- >>> :t dot sum (*) v v +-- dot sum (*) v v :: Array '[] Int +-- +-- >>> dot sum (*) v v +-- 14 +-- +-- dot sum (*) m v on a matrix and a vector is matrix-vector multiplication +-- Note that an `Array '[3] Int` is neither a row vector nor column vector. `dot` is not turning the vector into a matrix and then using matrix multiplication. +-- +-- >>> dot sum (*) v b +-- [9, 12, 15] +-- +-- >>> dot sum (*) b v +-- [14, 32] +dot :: + forall a b c d sa sb s' ss se. + ( HasShape sa, + HasShape sb, + HasShape (sa ++ sb), + se ~ TakeIndexes (sa ++ sb) '[Rank sa - 1, Rank sa], + HasShape se, + KnownNat (Minimum se), + KnownNat (Rank sa - 1), + KnownNat (Rank sa), + ss ~ '[Minimum se], + HasShape ss, + s' ~ DropIndexes (sa ++ sb) '[Rank sa - 1, Rank sa], + HasShape s' + ) => + (Array ss c -> d) -> + (a -> b -> c) -> + Array sa a -> + Array sb b -> + Array s' d +dot f g a b = contract f (Proxy :: Proxy '[Rank sa - 1, Rank sa]) (expand g a b) + +-- | select elements along every dimension +-- +-- >>> let s = slice (Proxy :: Proxy '[[0,1],[0,2],[1,2]]) a +-- >>> :t s +-- s :: Array '[2, 2, 2] Int +-- +-- >>> s +-- [[[2, 3], +-- [10, 11]], +-- [[14, 15], +-- [22, 23]]] +-- +-- >>> let s = squeeze $ slice (Proxy :: Proxy '[ '[0], '[0], '[0]]) a +-- >>> :t s +-- s :: Array '[] Int +-- +-- >>> s +-- 1 +slice :: + forall (pss :: [[Nat]]) s s' a. + ( HasShape s, + HasShape s', + KnownNatss pss, + KnownNat (Rank pss), + s' ~ Ranks pss + ) => + Proxy pss -> + Array s a -> + Array s' a +slice pss a = tabulate go + where + go s = index a (zipWith (!!) pss' s) + pss' = natValss pss + +-- | remove singleton dimensions +-- +-- >>> let a = [1..24] :: Array '[2,1,3,4,1] Int +-- >>> a +-- [[[[[1], +-- [2], +-- [3], +-- [4]], +-- [[5], +-- [6], +-- [7], +-- [8]], +-- [[9], +-- [10], +-- [11], +-- [12]]]], +-- [[[[13], +-- [14], +-- [15], +-- [16]], +-- [[17], +-- [18], +-- [19], +-- [20]], +-- [[21], +-- [22], +-- [23], +-- [24]]]]] +-- >>> squeeze a +-- [[[1, 2, 3, 4], +-- [5, 6, 7, 8], +-- [9, 10, 11, 12]], +-- [[13, 14, 15, 16], +-- [17, 18, 19, 20], +-- [21, 22, 23, 24]]] +-- +-- >>> squeeze ([1] :: Array '[1,1] Double) +-- 1.0 +squeeze :: + forall s t a. + (t ~ Squeeze s) => + Array s a -> + Array t a +squeeze (Array x) = Array x + +-- $scalar +-- Scalar specialisations + +-- | +-- +-- An /Array '[] a/ despite being a Scalar is never-the-less a one-element vector under the hood. Unification of representation is unexplored. +type Scalar a = Array ('[] :: [Nat]) a + +-- | unwrapping scalars is probably a performance bottleneck +-- +-- >>> let s = [3] :: Array ('[] :: [Nat]) Int +-- >>> fromScalar s +-- 3 +fromScalar :: (HasShape ('[] :: [Nat])) => Array ('[] :: [Nat]) a -> a +fromScalar a = index a ([] :: [Int]) + +-- | convert a number to a scalar +-- +-- >>> :t toScalar 2 +-- toScalar 2 :: Num a => Array '[] a +toScalar :: (HasShape ('[] :: [Nat])) => a -> Array ('[] :: [Nat]) a +toScalar a = fromList [a] + +-- $vector +-- Vector specialisations + +-- | +type Vector s a = Array '[s] a + +-- $matrix +-- Matrix specialisations + +-- | +type Matrix m n a = Array '[m, n] a + +instance + ( Multiplicative a, + P.Distributive a, + Subtractive a, + KnownNat m, + HasShape '[m, m] + ) => + Multiplicative (Matrix m m a) + where + + (*) = mmult + + one = ident + +-- | extract specialised to a matrix +-- +-- >>> row 1 m +-- [4, 5, 6, 7] +row :: forall m n a. (KnownNat m, KnownNat n, HasShape '[m, n]) => Int -> Matrix m n a -> Vector n a +row i (Array a) = Array $ V.slice (i * n) n a + where + n = fromIntegral $ natVal @n Proxy + +-- | row extraction checked at type level +-- +-- >>> safeRow (Proxy :: Proxy 1) m +-- [4, 5, 6, 7] +-- +-- >>> safeRow (Proxy :: Proxy 3) m +-- ... +-- ... index outside range +-- ... +safeRow :: forall m n a j. ('True ~ CheckIndex j m, KnownNat j, KnownNat m, KnownNat n, HasShape '[m, n]) => Proxy j -> Matrix m n a -> Vector n a +safeRow _j (Array a) = Array $ V.slice (j * n) n a + where + n = fromIntegral $ natVal @n Proxy + j = fromIntegral $ natVal @j Proxy + +-- | extract specialised to a matrix +-- +-- >>> col 1 m +-- [1, 5, 9] +col :: forall m n a. (KnownNat m, KnownNat n, HasShape '[m, n]) => Int -> Matrix m n a -> Vector n a +col i (Array a) = Array $ V.generate m (\x -> V.unsafeIndex a (i + x * n)) + where + m = fromIntegral $ natVal @m Proxy + n = fromIntegral $ natVal @n Proxy + +-- | column extraction checked at type level +-- +-- >>> safeCol (Proxy :: Proxy 1) m +-- [1, 5, 9] +-- +-- >>> safeCol (Proxy :: Proxy 4) m +-- ... +-- ... index outside range +-- ... +safeCol :: forall m n a j. ('True ~ CheckIndex j n, KnownNat j, KnownNat m, KnownNat n, HasShape '[m, n]) => Proxy j -> Matrix m n a -> Vector n a +safeCol _j (Array a) = Array $ V.generate m (\x -> V.unsafeIndex a (j + x * n)) + where + m = fromIntegral $ natVal @m Proxy + n = fromIntegral $ natVal @n Proxy + j = fromIntegral $ natVal @j Proxy + +-- | matrix multiplication +-- +-- This is dot sum (*) specialised to matrices +-- +-- >>> let a = [1, 2, 3, 4] :: Array '[2, 2] Int +-- >>> let b = [5, 6, 7, 8] :: Array '[2, 2] Int +-- >>> a +-- [[1, 2], +-- [3, 4]] +-- +-- >>> b +-- [[5, 6], +-- [7, 8]] +-- +-- >>> mmult a b +-- [[19, 22], +-- [43, 50]] +mmult :: + forall m n k a. + ( KnownNat k, + KnownNat m, + KnownNat n, + HasShape [m, n], + Ring a + ) => + Array [m, k] a -> + Array [k, n] a -> + Array [m, n] a +mmult (Array x) (Array y) = tabulate go + where + go [] = throw (NumHaskException "Needs two dimensions") + go [_] = throw (NumHaskException "Needs two dimensions") + go (i : j : _) = sum $ V.zipWith (*) (V.slice (fromIntegral i * k) k x) (V.generate k (\x' -> y V.! (fromIntegral j + x' * n))) + n = fromIntegral $ natVal @n Proxy + k = fromIntegral $ natVal @k Proxy +{-# INLINE mmult #-} diff --git a/src/NumHask/Array/HMatrix.hs b/src/NumHask/Array/HMatrix.hs new file mode 100644 --- /dev/null +++ b/src/NumHask/Array/HMatrix.hs @@ -0,0 +1,149 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE DeriveGeneric #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE NoStarIsType #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeApplications #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# OPTIONS_GHC -Wno-redundant-constraints #-} +{-# OPTIONS_GHC -fno-warn-incomplete-uni-patterns #-} +{-# OPTIONS_GHC -fno-warn-missing-signatures #-} + +module NumHask.Array.HMatrix + ( Array(..), + index, + mmult, + ) where + +import GHC.Exts (IsList (..)) +import GHC.Show (Show (..)) +import NumHask.Array.Shape +import NumHask.Prelude as P +import qualified Numeric.LinearAlgebra as H +import qualified Numeric.LinearAlgebra.Devel as H +import qualified Prelude + +newtype Array s a = Array {unArray :: H.Matrix a} + deriving (Show, NFData, Generic) + +-- | explicit rather than via Representable +index :: + forall s a. + ( HasShape s, + H.Element a, + H.Container H.Vector a + ) => + Array s a -> + [Int] -> + a +index (Array v) i = H.flatten v `H.atIndex` flatten s i + where + s = shapeVal (toShape @s) + +instance + ( Additive a, + HasShape s, + H.Container H.Vector a, + Num a + ) => + Additive (Array s a) + where + + (+) (Array x1) (Array x2) = Array $ H.add x1 x2 + + zero = Array $ H.konst zero (n, m) + where + s = shapeVal (toShape @s) + [n, m] = s + +instance + ( Multiplicative a, + HasShape s, + H.Container H.Vector a, + Num (H.Vector a), + Num a + ) => + Multiplicative (Array s a) + where + + (*) (Array x1) (Array x2) = Array $ H.liftMatrix2 (Prelude.*) x1 x2 + + one = Array $ H.konst one (n, m) + where + s = shapeVal (toShape @s) + [n, m] = s + +type instance Actor (Array s a) = a + +instance + ( HasShape s, + P.Distributive a, + CommutativeRing a, + Semiring a, + H.Container H.Vector a, + Num (H.Vector a), + Num a + ) => + Hilbert (Array s a) + where + (<.>) (Array a) (Array b) = H.sumElements $ H.liftMatrix2 (Prelude.*) a b + {-# INLINE (<.>) #-} + +instance + ( HasShape s, + Multiplicative a, + H.Container H.Vector a, + Num a + ) => + MultiplicativeAction (Array s a) + where + + (.*) (Array r) s = Array $ H.cmap (* s) r + {-# INLINE (.*) #-} + + (*.) s (Array r) = Array $ H.cmap (s *) r + {-# INLINE (*.) #-} + +-- | from flat list +instance + ( HasShape s, + Additive a, + H.Element a + ) => + IsList (Array s a) + where + + type Item (Array s a) = a + + fromList l = Array $ H.reshape n $ H.fromList $ take mn $ l ++ repeat zero + where + mn = P.product $ shapeVal (toShape @s) + s = shapeVal (toShape @s) + n = Prelude.last s + + toList (Array v) = H.toList $ H.flatten v + +-- | fast +mmult :: + forall m n k a. + ( KnownNat k, + KnownNat m, + KnownNat n, + HasShape [m, n], + Ring a, + H.Numeric a + ) => + Array [m, k] a -> + Array [k, n] a -> + Array [m, n] a +mmult (Array x) (Array y) = Array $ x H.<> y + +type instance Actor (Array s a) = a diff --git a/src/NumHask/Array/Shape.hs b/src/NumHask/Array/Shape.hs new file mode 100644 --- /dev/null +++ b/src/NumHask/Array/Shape.hs @@ -0,0 +1,399 @@ +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE NoStarIsType #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeApplications #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeInType #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE UndecidableInstances #-} +{-# OPTIONS_GHC -Wall #-} + +-- | Functions for manipulating shape. The module tends to supply equivalent functionality at type-level and value-level with functions of the same name (except for capitalization). +module NumHask.Array.Shape where + +import Data.List ((!!)) +import Data.Type.Bool +import GHC.TypeLits as L +import NumHask.Prelude as P hiding (Last, minimum) + +-- | The Shape type holds a [Nat] at type level and the equivalent [Int] at value level. +newtype Shape (s :: [Nat]) = Shape {shapeVal :: [Int]} deriving (Show) + +class HasShape s where + toShape :: Shape s + +instance HasShape '[] where + toShape = Shape [] + +instance (KnownNat n, HasShape s) => HasShape (n : s) where + toShape = Shape $ fromInteger (natVal (Proxy :: Proxy n)) : shapeVal (toShape :: Shape s) + +-- | Number of dimensions +rank :: [a] -> Int +rank = length +{-# INLINE rank #-} + +type family Rank (s :: [a]) :: Nat where + Rank '[] = 0 + Rank (_ : s) = Rank s + 1 + +-- | The shape of a list of element indexes +ranks :: [[a]] -> [Int] +ranks = fmap rank +{-# INLINE ranks #-} + +type family Ranks (s :: [[a]]) :: [Nat] where + Ranks '[] = '[] + Ranks (x : xs) = Rank x : Ranks xs + +-- | Number of elements +size :: [Int] -> Int +size [] = 1 +size [x] = x +size xs = P.product xs +{-# INLINE size #-} + +type family Size (s :: [Nat]) :: Nat where + Size '[] = 1 + Size (n : s) = n L.* Size s + +-- | convert from n-dim shape index to a flat index +-- +-- >>> flatten [2,3,4] [1,1,1] +-- 17 +-- +-- >>> flatten [] [1,1,1] +-- 0 +flatten :: [Int] -> [Int] -> Int +flatten [] _ = 0 +flatten _ [x'] = x' +flatten ns xs = sum $ zipWith (*) xs (drop 1 $ scanr (*) one ns) +{-# INLINE flatten #-} + +-- | convert from a flat index to a shape index +-- +-- >>> shapen [2,3,4] 17 +-- [1,1,1] +shapen :: [Int] -> Int -> [Int] +shapen [] _ = [] +shapen [_] x' = [x'] +shapen [_, y] x' = let (i, j) = divMod x' y in [i, j] +shapen ns x = + fst $ + foldr + ( \a (acc, r) -> + let (d, m) = divMod r a + in (m : acc, d) + ) + ([], x) + ns +{-# INLINE shapen #-} + +-- | /checkIndex i n/ checks if /i/ is a valid index of a list of length /n/ +checkIndex :: Int -> Int -> Bool +checkIndex i n = zero <= i && i + one <= n + +type family CheckIndex (i :: Nat) (n :: Nat) :: Bool where + CheckIndex i n = + If ((0 <=? i) && (i + 1 <=? n)) 'True (L.TypeError ('Text "index outside range")) + +-- | /checkIndexes is n/ check if /is/ are valid indexes of a list of length /n/ +checkIndexes :: [Int] -> Int -> Bool +checkIndexes is n = all (`checkIndex` n) is + +type family CheckIndexes (i :: [Nat]) (n :: Nat) :: Bool where + CheckIndexes '[] n = 'True + CheckIndexes (i : is) n = CheckIndex i n && CheckIndexes is n + +-- | dimension i is the i'th dimension of a Shape +dimension :: [Int] -> Int -> Int +dimension (s : _) 0 = s +dimension (_ : s) n = dimension s (n - 1) +dimension _ _ = throw (NumHaskException "dimension overflow") + +type family Dimension (s :: [Nat]) (i :: Nat) :: Nat where + Dimension (s : _) 0 = s + Dimension (_ : s) n = Dimension s (n - 1) + Dimension _ _ = L.TypeError ('Text "dimension overflow") + +-- | minimum value in a list +minimum :: [Int] -> Int +minimum [] = throw (NumHaskException "dimension underflow") +minimum [x] = x +minimum (x : xs) = P.min x (minimum xs) + +type family Minimum (s :: [Nat]) :: Nat where + Minimum '[] = L.TypeError ('Text "zero dimension") + Minimum '[x] = x + Minimum (x : xs) = If (x <=? Minimum xs) x (Minimum xs) + +type family Take (n :: Nat) (a :: [k]) :: [k] where + Take 0 _ = '[] + Take n (x : xs) = x : Take (n - 1) xs + +type family Drop (n :: Nat) (a :: [k]) :: [k] where + Drop 0 xs = xs + Drop n (_ : xs) = Drop (n - 1) xs + +type family Tail (a :: [k]) :: [k] where + Tail '[] = L.TypeError ('Text "No tail") + Tail (_ : xs) = xs + +type family Init (a :: [k]) :: [k] where + Init '[] = L.TypeError ('Text "No init") + Init '[_] = '[] + Init (x : xs) = x : Init xs + +type family Head (a :: [k]) :: k where + Head '[] = L.TypeError ('Text "No head") + Head (x : _) = x + +type family Last (a :: [k]) :: k where + Last '[] = L.TypeError ('Text "No last") + Last '[x] = x + Last (_ : xs) = Last xs + +type family (a :: [k]) ++ (b :: [k]) :: [k] where + '[] ++ b = b + (a : as) ++ b = a : (as ++ b) + +-- | drop the i'th dimension from a shape +-- +-- >>> dropIndex [2, 3, 4] 1 +-- [2,4] +dropIndex :: [Int] -> Int -> [Int] +dropIndex s i = take i s ++ drop (i + 1) s + +type DropIndex s i = Take i s ++ Drop (i + 1) s + +-- | /addIndex s i d/ adds a new dimension to shape /s/ at position /i/ +-- +-- >>> addIndex [2,4] 1 3 +-- [2,3,4] +addIndex :: [Int] -> Int -> Int -> [Int] +addIndex s i d = take i s ++ (d : drop i s) + +type AddIndex s i d = Take i s ++ (d : Drop i s) + +type Reverse (a :: [k]) = ReverseGo a '[] + +type family ReverseGo (a :: [k]) (b :: [k]) :: [k] where + ReverseGo '[] b = b + ReverseGo (a : as) b = ReverseGo as (a : b) + +-- | convert a list of position that references a final shape to one that references positions relative to an accumulator. Deletions are from the left and additions are from the right. +-- +-- deletions +-- +-- >>> posRelative [0,1] +-- [0,0] +-- +-- additions +-- +-- >>> reverse (posRelative (reverse [1,0])) +-- [0,0] +posRelative :: [Int] -> [Int] +posRelative as = reverse (go [] as) + where + go r [] = r + go r (x : xs) = go (x : r) ((\y -> bool (y - one) y (y < x)) <$> xs) + +type family PosRelative (s :: [Nat]) where + PosRelative s = PosRelativeGo s '[] + +type family PosRelativeGo (r :: [Nat]) (s :: [Nat]) where + PosRelativeGo '[] r = Reverse r + PosRelativeGo (x : xs) r = PosRelativeGo (DecMap x xs) (x : r) + +type family DecMap (x :: Nat) (ys :: [Nat]) :: [Nat] where + DecMap _ '[] = '[] + DecMap x (y : ys) = If (y + 1 <=? x) y (y - 1) : DecMap x ys + +-- | drop dimensions of a shape according to a list of positions (where position refers to the initial shape) +-- +-- >>> dropIndexes [2, 3, 4] [1, 0] +-- [4] +dropIndexes :: [Int] -> [Int] -> [Int] +dropIndexes s i = foldl' dropIndex s (posRelative i) + +type family DropIndexes (s :: [Nat]) (i :: [Nat]) where + DropIndexes s i = DropIndexesGo s (PosRelative i) + +type family DropIndexesGo (s :: [Nat]) (i :: [Nat]) where + DropIndexesGo s '[] = s + DropIndexesGo s (i : is) = DropIndexesGo (DropIndex s i) is + +-- | insert a list of dimensions according to position and dimension lists. Note that the list of positions references the final shape and not the initial shape. +-- +-- >>> addIndexes [4] [1,0] [3,2] +-- [2,3,4] +addIndexes :: () => [Int] -> [Int] -> [Int] -> [Int] +addIndexes as xs = addIndexesGo as (reverse (posRelative (reverse xs))) + where + addIndexesGo as' [] _ = as' + addIndexesGo as' (x : xs') (y : ys') = addIndexesGo (addIndex as' x y) xs' ys' + addIndexesGo _ _ _ = throw (NumHaskException "mismatched ranks") + +type family AddIndexes (as :: [Nat]) (xs :: [Nat]) (ys :: [Nat]) where + AddIndexes as xs ys = AddIndexesGo as (Reverse (PosRelative (Reverse xs))) ys + +type family AddIndexesGo (as :: [Nat]) (xs :: [Nat]) (ys :: [Nat]) where + AddIndexesGo as' '[] _ = as' + AddIndexesGo as' (x : xs') (y : ys') = AddIndexesGo (AddIndex as' x y) xs' ys' + AddIndexesGo _ _ _ = L.TypeError ('Text "mismatched ranks") + +-- | take list of dimensions according to position lists. +-- +-- >>> takeIndexes [2,3,4] [2,0] +-- [4,2] +takeIndexes :: [Int] -> [Int] -> [Int] +takeIndexes s i = (s !!) <$> i + +type family TakeIndexes (s :: [Nat]) (i :: [Nat]) where + TakeIndexes '[] _ = '[] + TakeIndexes _ '[] = '[] + TakeIndexes s (i : is) = + (s !! i) ': TakeIndexes s is + +type family (a :: [k]) !! (b :: Nat) :: k where + (!!) '[] i = L.TypeError ('Text "Index Underflow") + (!!) (x : _) 0 = x + (!!) (_ : xs) i = (!!) xs (i - 1) + +type family Enumerate (n :: Nat) where + Enumerate n = Reverse (EnumerateGo n) + +type family EnumerateGo (n :: Nat) where + EnumerateGo 0 = '[] + EnumerateGo n = (n - 1) : EnumerateGo (n - 1) + +-- | turn a list of included positions for a given rank into a list of excluded positions +-- +-- >>> exclude 3 [1,2] +-- [0] +exclude :: Int -> [Int] -> [Int] +exclude r = dropIndexes [0 .. (r - 1)] + +type family Exclude (r :: Nat) (i :: [Nat]) where + Exclude r i = DropIndexes (EnumerateGo r) i + +concatenate' :: Int -> [Int] -> [Int] -> [Int] +concatenate' i s0 s1 = take i s0 ++ (dimension s0 i + dimension s1 i : drop (i + 1) s0) + +type Concatenate i s0 s1 = Take i s0 ++ (Dimension s0 i + Dimension s1 i : Drop (i + 1) s0) + +type CheckConcatenate i s0 s1 s = + ( CheckIndex i (Rank s0) + && DropIndex s0 i == DropIndex s1 i + && Rank s0 == Rank s1 + ) + ~ 'True + +type CheckInsert d i s = + (CheckIndex d (Rank s) && CheckIndex i (Dimension s d)) ~ 'True + +type Insert d s = Take d s ++ (Dimension s d + 1 : Drop (d + 1) s) + +-- | /incAt d s/ increments the index at /d/ of shape /s/ by one. +incAt :: Int -> [Int] -> [Int] +incAt d s = take d s ++ (dimension s d + 1 : drop (d + 1) s) + +-- | /decAt d s/ decrements the index at /d/ of shape /s/ by one. +decAt :: Int -> [Int] -> [Int] +decAt d s = take d s ++ (dimension s d - 1 : drop (d + 1) s) + +-- /reorder' s i/ reorders the dimensions of shape /s/ according to a list of positions /i/ +-- +-- >>> reorder' [2,3,4] [2,0,1] +-- [4,2,3] +reorder' :: [Int] -> [Int] -> [Int] +reorder' [] _ = [] +reorder' _ [] = [] +reorder' s (d : ds) = dimension s d : reorder' s ds + +type family Reorder (s :: [Nat]) (ds :: [Nat]) :: [Nat] where + Reorder '[] _ = '[] + Reorder _ '[] = '[] + Reorder s (d : ds) = Dimension s d : Reorder s ds + +type family CheckReorder (ds :: [Nat]) (s :: [Nat]) where + CheckReorder ds s = + If + ( Rank ds == Rank s + && CheckIndexes ds (Rank s) + ) + 'True + (L.TypeError ('Text "bad dimensions")) + ~ 'True + +squeeze' :: (Eq a, Num a) => [a] -> [a] +squeeze' = filter (/=1) + +type family Squeeze (a :: [Nat]) where + Squeeze '[] = '[] + Squeeze a = Filter '[] a 1 + +type family Filter (r :: [Nat]) (xs :: [Nat]) (i :: Nat) where + Filter r '[] _ = Reverse r + Filter r (x : xs) i = Filter (If (x == i) r (x : r)) xs i + +-- unused but useful type-level functions + +type family Sort (xs :: [k]) :: [k] where + Sort '[] = '[] + Sort (x ': xs) = (Sort (SFilter 'FMin x xs) ++ '[x]) ++ Sort (SFilter 'FMax x xs) + +data Flag = FMin | FMax + +type family Cmp (a :: k) (b :: k) :: Ordering + +type family SFilter (f :: Flag) (p :: k) (xs :: [k]) :: [k] where + SFilter f p '[] = '[] + SFilter 'FMin p (x ': xs) = If (Cmp x p == 'LT) (x ': SFilter 'FMin p xs) (SFilter 'FMin p xs) + SFilter 'FMax p (x ': xs) = If (Cmp x p == 'GT || Cmp x p == 'EQ) (x ': SFilter 'FMax p xs) (SFilter 'FMax p xs) + +type family Zip lst lst' where + Zip lst lst' = ZipWith '(,) lst lst' -- Implemented as TF because #11375 + +type family ZipWith f lst lst' where + ZipWith f '[] lst = '[] + ZipWith f lst '[] = '[] + ZipWith f (l ': ls) (n ': ns) = f l n ': ZipWith f ls ns + +type family Fst a where + Fst '(a, _) = a + +type family Snd a where + Snd '(_, a) = a + +type family FMap f lst where + FMap f '[] = '[] + FMap f (l ': ls) = f l ': FMap f ls + +-- | Reflect a list of Nats +class KnownNats (ns :: [Nat]) where + natVals :: Proxy ns -> [Int] + +instance KnownNats '[] where + natVals _ = [] + +instance (KnownNat n, KnownNats ns) => KnownNats (n : ns) where + natVals _ = fromInteger (natVal (Proxy @n)) : natVals (Proxy @ns) + +-- | Reflect a list of list of Nats +class KnownNatss (ns :: [[Nat]]) where + natValss :: Proxy ns -> [[Int]] + +instance KnownNatss '[] where + natValss _ = [] + +instance (KnownNats n, KnownNatss ns) => KnownNatss (n : ns) where + natValss _ = natVals (Proxy @n) : natValss (Proxy @ns) diff --git a/src/NumHask/Shape.hs b/src/NumHask/Shape.hs deleted file mode 100644 --- a/src/NumHask/Shape.hs +++ /dev/null @@ -1,22 +0,0 @@ -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE UndecidableInstances #-} -{-# OPTIONS_GHC -Wall #-} - --- | numbers with a shape -module NumHask.Shape - ( HasShape(..) - -- * Representable - -- | Representable has most of what's needed to define numbers that have elements (aka scalars) and a fixed shape. - , Representable(..) - ) -where - -import Data.Functor.Rep - --- | Not everything that has a shape is representable. --- --- todo: Structure is a useful alternative concept/naming convention -class HasShape f where - type Shape f - shape :: f a -> Shape f diff --git a/test/test.hs b/test/test.hs --- a/test/test.hs +++ b/test/test.hs @@ -1,105 +1,66 @@ -{-# LANGUAGE CPP #-} -{-# LANGUAGE ConstraintKinds #-} {-# LANGUAGE DataKinds #-} -{-# LANGUAGE ExistentialQuantification #-} {-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE GADTs #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE NoImplicitPrelude #-} -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE PolyKinds #-} -{-# LANGUAGE RankNTypes #-} -{-# LANGUAGE RebindableSyntax #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeApplications #-} {-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE UndecidableInstances #-} {-# OPTIONS_GHC -Wall #-} -{-# OPTIONS_GHC -Wno-redundant-constraints #-} module Main where -import Data.Functor.Rep -import GHC.Exts (IsList(..)) -import NumHask.Array +import GHC.Exts (IsList (..)) +import qualified Hedgehog as H +import NumHask.Array.Fixed +import NumHask.Array.Shape import NumHask.Hedgehog import NumHask.Prelude as P -import Numeric.Dimensions as D import Test.DocTest -import qualified Hedgehog as H -import qualified NumHask.Hedgehog.Prop.Space as I import qualified Prelude -genAIntegral :: forall a m r. (H.MonadGen m, Dimensions r, Additive a, Bounded a, ToInteger a, FromInteger a) => m (Array [] (r :: [Nat]) a) +genAIntegral :: forall a m r. (HasShape r, H.MonadGen m, Additive a, Bounded a, ToInteger a, FromInteger a) => m (Array (r :: [Nat]) a) genAIntegral = fromList <$> replicateM (fromIntegral n) integral_ where - n = totalDim $ dims @Nat @r + n = product $ shapeVal $ toShape @r -genARational :: forall a m r. (H.MonadGen m, Dimensions r, Field a, Subtractive a, ToRatio a, FromRatio a) => m (Array [] (r :: [Nat]) a) +genARational :: forall a m r. (H.MonadGen m, HasShape r, Field a, Subtractive a, ToRatio a Integer, FromRatio a Integer) => m (Array (r :: [Nat]) a) genARational = fromList <$> replicateM (fromIntegral n) negUniform where - n = totalDim $ dims @Nat @r + n = product $ shapeVal $ toShape @r main :: IO () main = do - putStrLn ("Array DocTest turned on" :: Text) - doctest ["src/NumHask/Array.hs"] - putStrLn ("Example DocTest turned on" :: Text) - doctest ["src/NumHask/Array/Example.hs"] - bVInt <- assertProps "Vector Int 6" (Prelude.fromInteger 100) - (genAIntegral :: H.Gen (Vector [] 6 Int)) integralProps' - bMInt <- assertProps "Matrix [] '[3,4] Int" (Prelude.fromInteger 100) - (genAIntegral :: H.Gen (Array [] '[3,4] Int)) integralProps' - -- bVFloat <- assertProps "Vector Float 6" (Prelude.fromInteger 100) - -- (genARational :: H.Gen (Vector [] 6 Float)) (fieldProps' acc) - bMFloat <- assertProps "Array [] '[3,4] Float" (Prelude.fromInteger 100) - (genARational :: H.Gen (Array [] '[3,4] Float)) (fieldProps' acc) - unless (bVInt && bMInt && bMFloat) + putStrLn ("NumHask.Array.Fixed DocTest turned on" :: Text) + doctest ["src/NumHask/Array/Fixed.hs"] + putStrLn ("NumHask.Array.Shape DocTest turned on" :: Text) + doctest ["src/NumHask/Array/Shape.hs"] + putStrLn ("NumHask.Array.Dynamic DocTest turned on" :: Text) + doctest ["src/NumHask/Array/Dynamic.hs"] + bVInt <- + assertProps + "Vector Int 6" + (Prelude.fromInteger 100) + (genAIntegral :: H.Gen (Vector 6 Int)) + numhaskProps + bMInt <- + assertProps + "Matrix '[3,4] Int" + (Prelude.fromInteger 100) + (genAIntegral :: H.Gen (Array '[3, 4] Int)) + numhaskProps + unless + (bVInt && bMInt) exitFailure - where - acc = tabulate (const 1.0) -integralProps' - :: forall a. - ( Show a - , Eq a - , Distributive a - , Subtractive a - , Signed a - ) - => H.Gen a - -> [(H.PropertyName, H.Property)] -integralProps' g = mconcat $ - (\x -> x g) <$> - [ isAdditive - , isSubtractive - , isMultiplicative - , \x -> [("distributive", isDistributive zero (+) (*) x)] - , \x -> [("signed", NumHask.Hedgehog.isSigned x)] - ] - --- | field laws -fieldProps' - :: forall a. - ( Show a - , Epsilon a - , Lattice a - , LowerBoundedField a - , BoundedJoinSemiLattice a - , BoundedMeetSemiLattice a - , Signed a - ) - => a - -> H.Gen a - -> [(H.PropertyName, H.Property)] -fieldProps' acc g = mconcat $ - (\x -> x g) <$> - [ I.isAdditive acc - , \x -> [("subtractive", I.isSubtractive acc x)] - , I.isMultiplicative acc - , \x -> [("distributive", I.isDistributiveTimesPlus one x)] - , \x -> [("divisive", I.isDivisive one x)] - , \x -> [("signed", I.isSigned one x)] - ] +numhaskProps :: + forall a. + ( Show a, + Eq a, + Subtractive a + ) => + H.Gen a -> + [(H.PropertyName, H.Property)] +numhaskProps g = + mconcat $ + (\x -> x g) + <$> [ isAdditive, + isSubtractive + ]