diff --git a/Data/VectorSpace/Free.hs b/Data/VectorSpace/Free.hs --- a/Data/VectorSpace/Free.hs +++ b/Data/VectorSpace/Free.hs @@ -7,30 +7,38 @@ -- Stability : experimental -- Portability : portable -- -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE CPP #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE ConstrainedClassMethods #-} module Data.VectorSpace.Free ( -- * Supported types + -- ** Fixed low dimension -- | These come from the package. V0 , V1 , V2 , V3 , V4 + -- ** Arbitrary dimension + , Sequence, FinSuppSeq -- * The vector-space type classes -- ** General -- | These come from the package. , AffineSpace(..), AdditiveGroup(..) , VectorSpace(..), InnerSpace(..), HasBasis(..) + -- ** Small + , OneDimensional(..) -- ** Free , FiniteFreeSpace(..) ) where import Data.AffineSpace import Data.VectorSpace +import Data.VectorSpace.Free.FiniteSupportedSequence (FinSuppSeq) +import Data.VectorSpace.Free.Sequence (Sequence) import Data.Basis import Data.MemoTrie @@ -43,10 +51,18 @@ import Linear.V4 import qualified Linear.Affine as LA -import Control.Lens +import Control.Lens ((^.), FoldableWithIndex, ifoldr) +import qualified Data.Foldable as Foldable + import qualified Data.Vector.Unboxed as UArr +import qualified Data.Vector.Generic.Mutable as MArr +import Data.Ratio +import Foreign.C.Types (CFloat, CDouble) + +import GHC.Exts (IsList(..)) + vDecomp :: FoldableWithIndex (L.E v) v => v s -> [(L.E v, s)] vDecomp = ifoldr (\b s l -> (b,s):l) [] @@ -116,9 +132,44 @@ enumerate (V4T (V4 x y z w)) = [(L.ex, x), (L.ey, y), (L.ez, z), (L.ew, w)] +infixr 7 ^/^, ^/! +class (VectorSpace v, Fractional (Scalar v)) => OneDimensional v where + -- | Compare the (directed) length of two vectors. + (^/^) :: v -> v -> Maybe (Scalar v) + -- | Unsafe version of '^/^'. + (^/!) :: v -> v -> Scalar v + v^/!w = case v^/^w of + Just μ -> μ + Nothing -> 1/0 +instance OneDimensional Float where + _^/^0 = Nothing + x^/^y = Just $ x/y + x^/!y = x/y +instance OneDimensional Double where + _^/^0 = Nothing + x^/^y = Just $ x/y + x^/!y = x/y +instance OneDimensional CFloat where + _^/^0 = Nothing + x^/^y = Just $ x/y + x^/!y = x/y +instance OneDimensional CDouble where + _^/^0 = Nothing + x^/^y = Just $ x/y + x^/!y = x/y +instance Integral i => OneDimensional (Ratio i) where + _^/^0 = Nothing + x^/^y = Just $ x/y + x^/!y = x/y +instance (Eq r, Fractional r) => OneDimensional (V1 r) where + _^/^V1 0 = Nothing + V1 x^/^V1 y = Just $ x/y + V1 x^/!V1 y = x/y + + class (VectorSpace v, Num (Scalar v)) => FiniteFreeSpace v where {-# MINIMAL freeDimension, toFullUnboxVect, unsafeFromFullUnboxVect #-} freeDimension :: Functor p => p v -> Int @@ -180,3 +231,7 @@ (UArr.unsafeIndex v 1) (UArr.unsafeIndex v 2) (UArr.unsafeIndex v 3) + + + + diff --git a/Data/VectorSpace/Free/FiniteSupportedSequence.hs b/Data/VectorSpace/Free/FiniteSupportedSequence.hs new file mode 100644 --- /dev/null +++ b/Data/VectorSpace/Free/FiniteSupportedSequence.hs @@ -0,0 +1,108 @@ +-- | +-- Module : Data.VectorSpace.Free.FiniteSupportedSequence +-- Copyright : (c) Justus Sagemüller 2016 +-- License : GPL v3 +-- +-- Maintainer : (@) sagemueller $ geo.uni-koeln.de +-- Stability : experimental +-- Portability : portable +-- +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE ConstrainedClassMethods #-} + +module Data.VectorSpace.Free.FiniteSupportedSequence ( + FinSuppSeq (..) + ) where + +import Data.AffineSpace +import Data.VectorSpace +import Data.Basis + +import qualified Data.Foldable as Foldable + +import qualified Data.Vector.Unboxed as UArr +import qualified Data.Vector.Generic.Mutable as MArr + +import GHC.Exts (IsList(..)) + + + + +-- | The space of finitely-supported sequences is an /infinite/-dimensional space. +-- An vector of length /l/ is here understood as an infinite sequence that begins +-- with /l/ nonzero values, and continues with infinite zeroes. +-- +-- You may also consider this as the type that languages like Octave/Matlab +-- (as well as Haskell's library) +-- approximate with their “vectors”, with one important difference: there is +-- no such thing as a dimensional-mismatch error, since we consider all these vectors +-- as elements of the same infinite-dimensional space. Adding two different-size +-- vectors will simply zero-pad the shorter, and unlike in Matlab this behaviour extends +-- consequently to matrix multiplication etc. (defined in +-- ) +-- +-- Of course it /can/ make sense to constrain the dimension, but for this the +-- type system should be used, not runtime checks. +-- +-- (This is the same +-- behaviour that the library +-- gives to the standard list and vector types, but the problem there is that it +-- can't use unboxed arrays as these are not functors, but unboxing is crucial for +-- performance.) +newtype FinSuppSeq n = FinSuppSeq { getFiniteSeq :: UArr.Vector n } + + +{-# INLINE liftU2FSS #-} +liftU2FSS :: UArr.Unbox n => (n -> n -> n) -> FinSuppSeq n -> FinSuppSeq n -> FinSuppSeq n +-- Adapted from: +-- http://hackage.haskell.org/package/linear-1.20.5/docs/src/Linear.Vector.html#line-200 +liftU2FSS f (FinSuppSeq u) (FinSuppSeq v) = FinSuppSeq $ case compare lu lv of + LT | lu == 0 -> v + | otherwise -> UArr.modify + (\ w -> Foldable.forM_ [0..lu-1] $ + \i -> MArr.unsafeWrite w i $ f (UArr.unsafeIndex u i) + (UArr.unsafeIndex v i)) v + EQ -> UArr.zipWith f u v + GT | lv == 0 -> u + | otherwise -> UArr.modify + (\ w -> Foldable.forM_ [0..lv-1] $ + \i -> MArr.unsafeWrite w i $ f (UArr.unsafeIndex u i) + (UArr.unsafeIndex v i)) u + where lu = UArr.length u + lv = UArr.length v + + +instance (Num n, UArr.Unbox n) => AffineSpace (FinSuppSeq n) where + type Diff (FinSuppSeq n) = FinSuppSeq n + (.-.) = (^-^) + (.+^) = (^+^) + +instance (Num n, UArr.Unbox n) => AdditiveGroup (FinSuppSeq n) where + zeroV = FinSuppSeq $ UArr.empty + (^+^) = liftU2FSS (+) + negateV (FinSuppSeq v) = FinSuppSeq $ UArr.map negate v + +instance (Num n, UArr.Unbox n) => VectorSpace (FinSuppSeq n) where + type Scalar (FinSuppSeq n) = n + μ*^FinSuppSeq v = FinSuppSeq $ UArr.map (μ*) v + +instance (Num n, AdditiveGroup n, UArr.Unbox n) => InnerSpace (FinSuppSeq n) where + FinSuppSeq v<.>FinSuppSeq w = UArr.sum (UArr.zipWith (*) v w) + +instance (Num n, UArr.Unbox n) => HasBasis (FinSuppSeq n) where + type Basis (FinSuppSeq n) = Int + basisValue i = FinSuppSeq $ UArr.replicate i 0 `UArr.snoc` 1 + decompose = zip [0..] . toList + decompose' (FinSuppSeq v) i = maybe 0 id $ v UArr.!? i + +instance UArr.Unbox n => IsList (FinSuppSeq n) where + type Item (FinSuppSeq n) = n + fromListN l = FinSuppSeq . fromListN l + fromList = FinSuppSeq . fromList + toList = toList . getFiniteSeq + +instance (UArr.Unbox n, Show n) => Show (FinSuppSeq n) where + show = ("fromList "++) . show . toList diff --git a/Data/VectorSpace/Free/Sequence.hs b/Data/VectorSpace/Free/Sequence.hs new file mode 100644 --- /dev/null +++ b/Data/VectorSpace/Free/Sequence.hs @@ -0,0 +1,228 @@ +-- | +-- Module : Data.VectorSpace.Free.Sequence +-- Copyright : (c) Justus Sagemüller 2016 +-- License : GPL v3 +-- +-- Maintainer : (@) sagemueller $ geo.uni-koeln.de +-- Stability : experimental +-- Portability : portable +-- +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE ConstrainedClassMethods #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE DeriveFunctor #-} + +module Data.VectorSpace.Free.Sequence ( + Sequence, BoxSequence, GSequence (..), minimumChunkSize + ) where + +import Data.VectorSpace.Free.FiniteSupportedSequence + +import Data.AffineSpace +import Data.VectorSpace +import Data.Basis + +import qualified Data.Foldable as Foldable + +import qualified Data.Vector as Arr +import qualified Data.Vector.Unboxed as UArr +import qualified Data.Vector.Generic as GArr +import qualified Data.Vector.Generic.Mutable as MArr + +import GHC.Exts (IsList(..)) + + + + +-- | The space of possibly-infinite sequences, isomorphic to the space of all lists +-- but implemented more efficiently (with exponentially-growing chunks of unboxed data, +-- so the overhead is amortised). In other words, this is a type of spine-lazy +-- but element-strict arrays. +-- +-- This space is dual to 'FinSuppSeq', which is completely strict. +data GSequence array n + = Sequence { + sequenceHeads :: !(array n) + -- ^ Length must be at most 'minimumChunkSize' in the outer constructor and + -- double in each deeper layer. (Specification subject to future change!) + -- If the length at depth 𝑑 is less than 2^𝑑, the remaining entries are + -- treated as zeroes. + , sequenceRemain :: GSequence array n + } + | SoloChunk { + chunkOffset :: !Int + , soloChunk :: !(array n) + -- ^ Length plus offset must be at most 'minimumChunkSize' if this is + -- the outer constructor and can double in each deeper layer. + } + + +type Sequence = GSequence UArr.Vector +type BoxSequence = GSequence Arr.Vector + +deriving instance Functor BoxSequence + +reboxSequence :: (GArr.Vector a n, GArr.Vector b n) => GSequence a n -> GSequence b n +reboxSequence (SoloChunk o c) = SoloChunk o $ Arr.convert c +reboxSequence (Sequence h r) = Sequence (Arr.convert h) $ reboxSequence r + +mapSequence :: (UArr.Unbox n, UArr.Unbox m) => (n -> m) -> Sequence n -> Sequence m +mapSequence f (SoloChunk i₀ chunk) = SoloChunk i₀ (UArr.map f chunk) +mapSequence f (Sequence hd rm) = Sequence (UArr.map f hd) (mapSequence f rm) + +{-# INLINE liftU2Seq #-} +liftU2Seq :: UArr.Unbox n => n -> (n -> n -> n) -> Sequence n -> Sequence n -> Sequence n +liftU2Seq defaultVal f (Sequence hu ru) (Sequence hv rv) + = (`Sequence`liftU2Seq defaultVal f ru rv) $ case compare lu lv of +-- Adapted from: +-- http://hackage.haskell.org/package/linear-1.20.5/docs/src/Linear.Vector.html#line-200 + LT | lu == 0 -> hv + | otherwise -> UArr.modify + (\ w -> Foldable.forM_ [0..lu-1] $ + \i -> MArr.unsafeWrite w i $ f (UArr.unsafeIndex hu i) + (UArr.unsafeIndex hv i)) hv + EQ -> UArr.zipWith f hu hv + GT | lv == 0 -> hu + | otherwise -> UArr.modify + (\ w -> Foldable.forM_ [0..lv-1] $ + \i -> MArr.unsafeWrite w i $ f (UArr.unsafeIndex hu i) + (UArr.unsafeIndex hv i)) hu + where lu = UArr.length hu + lv = UArr.length hv +liftU2Seq defaultVal f (Sequence hu ru) (SoloChunk ov cv) + | lu == 0, lv == 0 = Sequence UArr.empty ru + | lu == 0 = Sequence (UArr.replicate ov defaultVal UArr.++ cv) ru + | lv == 0 = Sequence hu ru + | lu >= lv+ov = Sequence (UArr.modify + (\w -> Foldable.forM_ [0..lv-1] $ + \i -> MArr.unsafeWrite w (i+ov) + $ f (UArr.unsafeIndex hu (i+ov)) + (UArr.unsafeIndex cv i)) hu) + ru + | ov == 0 = Sequence (UArr.modify + (\w -> Foldable.forM_ [0..lu-1] $ + \i -> MArr.unsafeWrite w i + $ f (UArr.unsafeIndex hu i) + (UArr.unsafeIndex cv i)) cv) + ru + | otherwise = Sequence (UArr.take ov hu UArr.++ UArr.modify + (\w -> Foldable.forM_ [ov..lu-1] $ + \i -> MArr.unsafeWrite w i + $ f (UArr.unsafeIndex hu i) + (UArr.unsafeIndex cv (i-ov))) cv) + ru + where lu = UArr.length hu + lv = UArr.length cv +liftU2Seq defaultVal f (SoloChunk ou cu) (Sequence hv rv) + | lu == 0, lv == 0 = Sequence UArr.empty rv + | lu == 0 = Sequence hv rv + | lv == 0 = Sequence (UArr.replicate ou defaultVal UArr.++ cu) rv + | lu+ou <= lv = Sequence (UArr.modify + (\w -> Foldable.forM_ [0..lu-1] $ + \i -> MArr.unsafeWrite w (i+ou) + $ f (UArr.unsafeIndex cu i) + (UArr.unsafeIndex hv (i+ou))) hv) + rv + | ou == 0 = Sequence (UArr.modify + (\w -> Foldable.forM_ [0..lv-1] $ + \i -> MArr.unsafeWrite w i + $ f (UArr.unsafeIndex cu i) + (UArr.unsafeIndex hv i)) cu) + rv + | otherwise = Sequence (UArr.take ou hv UArr.++ UArr.modify + (\w -> Foldable.forM_ [ou..lv-1] $ + \i -> MArr.unsafeWrite w i + $ f (UArr.unsafeIndex cu (i-ou)) + (UArr.unsafeIndex hv i)) cu) + rv + where lu = UArr.length cu + lv = UArr.length hv +liftU2Seq _ f (SoloChunk ou cu) (SoloChunk ov cv) + | lu == 0 = SoloChunk ov cv + | lv == 0 = SoloChunk ou cu + | ou >= ov, ou+lu <= ov+lv + = SoloChunk ov $ UArr.modify + (\ w -> Foldable.forM_ [0..lu-1] $ + \i -> MArr.unsafeWrite w i $ f (UArr.unsafeIndex cu i) + (UArr.unsafeIndex cv (i+δo))) cv + | ou <= ov, ou+lu >= ov+lv + = SoloChunk ou $ UArr.modify + (\ w -> Foldable.forM_ [0..lv-1] $ + \i -> MArr.unsafeWrite w i $ f (UArr.unsafeIndex cu (i-δo)) + (UArr.unsafeIndex cv i)) cu + | ou >= ov + = SoloChunk ov $ UArr.take δo cv UArr.++ UArr.modify + (\ w -> Foldable.forM_ [δo..lv-1] $ + \i -> MArr.unsafeWrite w i $ f (UArr.unsafeIndex cu (i-δo)) + (UArr.unsafeIndex cv i)) cu + | ou <= ov + = SoloChunk ou $ UArr.take (-δo) cu UArr.++ UArr.modify + (\ w -> Foldable.forM_ [-δo..lu-1] $ + \i -> MArr.unsafeWrite w i $ f (UArr.unsafeIndex cu i) + (UArr.unsafeIndex cv (i+δo))) cv + where lu = UArr.length cu + lv = UArr.length cv + δo = ou-ov + + +instance (Num n, UArr.Unbox n) => AffineSpace (Sequence n) where + type Diff (Sequence n) = Sequence n + (.-.) = (^-^) + (.+^) = (^+^) + +instance (Num n, UArr.Unbox n) => AdditiveGroup (Sequence n) where + zeroV = SoloChunk 0 UArr.empty + (^+^) = liftU2Seq 1 (+) + negateV = mapSequence negate + +instance (Num n, UArr.Unbox n) => VectorSpace (Sequence n) where + type Scalar (Sequence n) = n + μ*^v = mapSequence (μ*) v + +instance (Num n, UArr.Unbox n) => HasBasis (Sequence n) where + type Basis (Sequence n) = Int + basisValue = bv minimumChunkSize + where bv chunkSize i + | i IsList (Sequence n) where + type Item (Sequence n) = n + fromListN = fln minimumChunkSize + where fln chunkSize l ns + | l>chunkSize = let (h,r) = splitAt chunkSize ns + in Sequence (UArr.fromList h) + $ fln (chunkSize*2) (l-chunkSize) r + | otherwise = SoloChunk 0 $ UArr.fromList ns + fromList = fln minimumChunkSize + where fln chunkSize ns + | null r = SoloChunk 0 $ UArr.fromList ns + | otherwise = Sequence (UArr.fromList h) + $ fln (chunkSize*2) r + where (h,r) = splitAt chunkSize ns + toList = tl minimumChunkSize + where tl _ (SoloChunk o c) = replicate o 0 ++ toList c + tl chunkSize (Sequence h r) + = toList h ++ replicate (chunkSize - UArr.length h) 0 + ++ tl (chunkSize*2) r + +instance (UArr.Unbox n, Show n, Num n) => Show (Sequence n) where + show = ("fromList "++) . show . toList + + + +minimumChunkSize :: Int +minimumChunkSize = 1 diff --git a/free-vector-spaces.cabal b/free-vector-spaces.cabal --- a/free-vector-spaces.cabal +++ b/free-vector-spaces.cabal @@ -2,7 +2,7 @@ -- documentation, see http://haskell.org/cabal/users-guide/ name: free-vector-spaces -version: 0.1.1.0 +version: 0.1.2.0 synopsis: Instantiate the classes from the vector-space package with types from linear description: The package offers efficient vector types — where vector means /element of a free vector space/, i.e. @@ -48,6 +48,8 @@ library exposed-modules: Data.VectorSpace.Free + , Data.VectorSpace.Free.FiniteSupportedSequence + , Data.VectorSpace.Free.Sequence -- other-modules: -- other-extensions: build-depends: base >=4.6 && <5.1,